Dual Solutions and Stability Analysis of Cu-H 2 O-Casson Nanofluid Convection past a Heated Stretching/Shrinking Slippery Sheet in a Porous Medium

In this study, we examined the impact of Cu-H 2 O nanoparticles on two-dimensional Casson nano ﬂ uid ﬂ ows past permeable stretching/shrinking sheet embedded in a Darcy-Forchheimer porous medium in the presence of slipperiness of surface, suction/injection, viscous dissipation, and convective heating. Using some realistic assumptions and appropriate similarity transformations, the governing nonlinear partial di ﬀ erential equations were formulated and transformed into a system of nonlinear ordinary di ﬀ erential equations and then numerically solved by using the shooting technique. Numerical results are displayed for dimensionless ﬂ uid velocity and temperature pro ﬁ les, skin friction, and the local Nusselt number. The impacts of di ﬀ erent governing physical parameters on these quantities are presented and discussed using graphs, tables, and a chart. For the speci ﬁ c range of shrinking sheet, the result shows that dual solutions exist, and temporal stability analysis is performed by introducing small disturbances to determine the stable solutions. It is detected that the upper branch solution is hydrodynamically stable and substantially realistic; however, the lower branch solution is unstable and physically unachievable. The ﬂ uid ﬂ ow stability is obtained by enhancing the suction, surface slipperiness, and viscous dissipation parameters. However, augmenting the values of the Casson factor, Cu-H 2 O nanoparticle volume fraction, porous medium, porous medium inertia, and convective heating parameters increases the blow-up stability of the ﬂ uid ﬂ ow. The rate of heat transfer enhances with the increment in the Casson factor, porous medium, porous medium inertia, suction, velocity ratio, nanoparticle volume fraction, and convective heating parameters, whereas it reduces as the slipperiness of the surface and viscous dissipation parameters rise. Increment of Cu-H 2 O nanoparticle volume fraction into the Casson ﬂ uid boosts the heat transfer enhancement rate higher for the shrinking sheet surface.


Introduction
For a few decades, researchers have given due attention to the study of non-Newtonian fluids due to their significant features in the fields of industrial processes and technological sciences, such as polymer engineering, certain separation processes, manufacturing of papers and foods, and petroleum drilling, according to Bhattacharyya [1]. For instance, drilling muds, synthetic lubricants, clay coating, biological fluids like blood, certain oils, paints, and sugar solutions are common cases of non-Newtonian fluid types. The fundamental equations of Navier-Stokes cannot momentarily define the characteristics of the non-Newtonian fluid flow field due to the complexity in the mathematical expression of the flow problem. For non-Newtonian fluid, several models are defined based on rheological qualities, such as Bulky, Casson, Eyring-Powell, Seely, Oldroyd-B, Maxwell, Oldroyd-A, Carreau, Jeffrey, and Burger. From these models, the Casson model [2] is the most important model for the suspension and blood properties in our daily life. Magnetohydrodynamic Casson fluid under the aligned magnetic field with inconstant thickness was investigated by Saravana et al. [3]. Zhang et al. [4] studied the heat transport characteristic flow of the Casson fluid with electroosmosis forces. Recently, heat transfer on the Casson nanofluid overflow past a Riga plate (one of the external agents used to limit the friction force and control the flow of fluid) was investigated by Upreti et al. [5].
For its frequent applications in industrial and engineering, the theme of flow nanofluids has attracted due attention. Tawari and Das [6] did an analysis on nanofluids. Gupta et al. [7] discovered the astonishing thermal properties of nanofluids and improved the causes of nanofluids' thermal conductivity. In many industrial processes, heating and cooling of fluids are fundamental demands such as power manufacturing and delivery. Frequently, there is a need to enhance the process of cooling of high-energy equipment. [8][9][10] worked on how to improve the thermal conductivity of fluids, so that the heat transfer of the fluid enhances. Eastman et al. [11] investigated enhancements in thermal conductivity of fluid using Cu nanofluids, where just a 0.3% volume fraction of 10 nm diameter Cu nanoparticles led to an increment of up to 40% in the ethylene glycol thermal conductivity. Moreover, the convection heat transfer of nanofluids was also investigated by different researchers, and based on that, significant improvement was reported in the heat transfer rate. Zubair et al. [12] did an analysis on the magnetohydrodynamic Casson nanofluid flow with entropy generation in the presence of viscous dissipation. Saeed et al. [13] demonstrated threedimensional Casson's nanofluid flow past a rotating inclined disk with the influence of heat absorption/production. Khan et al. [14] studied the impact of an induced magnetic field on mixed convective stagnation flow of TiO 2 -Cu-H 2 O hybrid nanofluid towards a stretchable sheet surface. Rizwana et al. [15] concentrated on the non-Newtonian fluid flow at an oblique stagnation point of Cu-H 2 O nanofluid flow, concluding that the Casson fluid parameter makes the fluid velocity faster and decreases the boundary layer while the temperature profile drops with non-Newtonian parameter β and nanoparticle volume fraction ϕ, and the system heats up by the impact of viscous dissipation. More discussion on heat transfer enhancement is detailed by Lund et al. [16] and Hussanan et al. [17].
Due to its enormous applications in areas of industry and engineering, according to Hayat et al. [18], the boundary layer flows of non-Newtonian fluids and the heat transfer properties due to a stretching sheet surface have attracted researchers' attention. Crane [19] investigated the viscous fluid flows because of linearly stretching surface problems. Tamoor et al. [20] studied the MHD flow of Casson's fluid past a stretching cylinder. In contrast to stretching sheets, Wang [21] considered shrinking sheets, which exhibit quite different properties from that of the stretching sheet flow. In the shrinking surface case, because of the free-range fluid flow happening in the boundary layer, no possible solution is obtained. Consequently, the addition of a sufficient amount of wall mass suction (which controls the vorticity produced because of the shrinking of the sheet within the boundary layer) by Miklavcic and Wang [22] or by adding stagnation flow by Wang [23] may guarantee the existence of a solution that is not unique. This leads to the existence of a nonunique solution for the system of governing differential equations. For the non-unique solutions obtained, we use a mathematical technique (called temporal stability analysis) that is used for testing temporal stability. Experimentally, the lower branch solution which is a part of the solution of the system differential equations cannot be produced and, hence, should be analyzed. Nazar et al. [24] examined a nanofluid stagnation-point flow over a shrinking sheet. Layek et al. [25] presented boundary layer stagnation-point flow of non-Newtonian fluids past a shrinking/stretching sheet and showed that as the value of the velocity ratio (shrinking/ stretching) parameter increases, the velocity and thermal boundary layer thicknesses decrease. Again, Mandal and Layek [26] investigated the unsteady magnetohydrodynamic (MHD) mixed convective Casson fluid flow over a flat surface in the presence of slip conditions and blowing/suction and obtained significant results. More flow stability analysis and testing for the existence of a dual solution were detailed in the literature [27][28][29][30].
Studies of fluid flow and heat transfer in a porous medium have attracted the attention of researchers towards fluid flow in porous media for the past several decades because of the wide range of applications, such as water developments in geothermal supplies and drying technologies. For simulating porous media there are the Darcy, non-Darcy, and nonequilibrium models. For low-velocity flow with weak porosity conditions, a pioneering semiempirical equation was developed by Darcy. Accordingly, Forchheimer [31] prophesied a modified equation known as the Darcy-Forchheimer equation by introducing quadratic terms in the governing momentum equation of the fluid flow. Harris et al. [32] studied mixed convection boundary-layer flow over a vertical surface through a porous medium using the Brinkman model with slip. The Brinkman-extended Darcy model was used by Dogonchi et al. [33] to explore the natural convection of a Cu-H 2 O nanofluid across a porous medium and obtain the Darcy number which has a direct relationship with the intensity of the convective flow across the medium. Shaw et al. [34] studied the application to brain dynamics of the impact of entropy generation on the Darcy-Forchheimer flow of MnFe 2 O 4 -Casson/water nanofluid because of a rotating disk. In their recent study, Upreti et al. [35] discussed the hydromagnetic stagnation point magnetite ferrofluid flow over a convectively heated shrinking/stretching permeable sheet surface which is exposed to injection/suction in a Darcy-Forchheimer porous medium. Joshi et al. [36] investigated the Darcy-Forchheimer flow model in a threedimensional case and the heat transfer phenomenon of H 2 O-CNT nanofluid on a two-way stretchable surface, and it was observed that with an enhancement in the Eckert numbers along the two directions, for temperature curves, two patterns were obtained, the initial temperature outlines rose, and after that, they diminished. Moreover, Tadesse et al. [37] studied the overall influence of blowing/suction A thorough observation of the aforementioned abovecited reviews and works reveals that no scientific analysis has been done to study the mutual effects of the slipperiness of the surface, nanoparticle volume fraction, thermal and viscous dissipation, porous media, and stretching/ shrinking surface on the rate of heat transfer enhancement in a flow into the boundary layer of the Cu-H 2 O-Casson nanofluids. Filling the gap in the work of Upreti et al. [35] (in which they consider MHD Newtonian flow), this paper considers non-Newtonian flow subjected to a slippery surface. Since in this type of fluid flow problem, dual solutions possibly exist due to stretching/shrinking surface, the main goal of this paper is to do stability analysis to determine physically reliable and stable solutions and enhancement of heat transfer as a result of adding the Cu-H 2 O nanoparticle to non-Newtonian Casson fluids for the coupled transformed ordinary differential equations. Moreover, the numerical results of the coefficient of skin friction, rates of heat transfer, and velocity and temperature profiles are demonstrated and discussed both graphically and quantitatively. Thermally driven boundary layer flow of the Casson nanofluid in a porous medium has several applications in chemical and mechanical engineering, e.g., food processing and storage, geophysical systems, electrochemistry, fibrous insulation, metallurgy, the design of pebble bed nuclear reactors, underground disposal of nuclear or nonnuclear waste, and microelectronics cooling. Casson's nanofluid is more helpful for cooling and friction-reducing agents compared to Newtonian-type nanofluid flow.

Mathematical Description of the Problem
Consider the two-dimensional, viscous, homogeneous, steady, isotropic, laminar, and incompressible boundary layer flow of a Cu-H 2 O-Casson nanofluid over a linearly stretching/shrinking sheet within a Darcy-Forchheimer porous medium. The flow is subjected to suction/injection of constant velocity V 0 perpendicular to the sheet surface. The shrinking/stretching sheet moves with velocity U w = ax, and the flow is subjected to slip condition. The surface below the shrinking/stretching sheet surface is heated convectively by a hot fluid having an initial temperature T f that gives a heat transfer coefficient h f , and the ambient temperature of the Casson nanofluid is taken as T ∞ . Geometrical orientation for the flow considers the Cartesian coordinate system where the x-axis is taken as the shrinking/stretching sheet surface, and the y-axis is perpendicular to the sheet surface so that the flow is confined to the half plane y > 0, as the physical model shown in Figure 1.
The rheological equation of an isotropic and incompressible flow of a Casson fluid can be written based on Nakamura and Sawada [38] and Animasaun [39] as where τ ij represents components of the stress tensor, μ B is the plastic dynamic viscosity of the non-Newtonian fluid, p y ≡ μ B ffiffiffiffiffiffiffi 2π c p /β is the yield stress of the fluid, β is the non-

Computational and Mathematical Methods
Newtonian Casson parameter, π = e ij e ij is the ði, jÞ th deformation rate component (product of the rate of strain tensor with itself), e ij = 1/2½ð∂u i /∂x j Þ/ð∂u j /∂x i Þ is the rate of strain tensor, and π c is a critical value of π, which is based on the non-Newtonian model. Some non-Newtonian fluids (for instance rheopectic fluids) require a gradually enhancing shear stress to get a constant strain rate. For the Casson fluid flow under consideration, π > π c , and the dynamic viscosity is defined as μ f = μ B + p y / ffiffiffiffiffi ffi 2π p , and hence, the kinematic viscosity becomes Adopting the Darcy-Forchheimer flow model, Tiwari-Das [6] convective transport model equations, for this investigation, the governing equations are formulated from the balance of continuity, linear momentum, and energy past a permeable shrinking/stretching sheet surface, with respect to a Cartesian coordinate; x-y system is given by Tshivhi and Makinde [28] and Tadesse et al. [37] as with the subjected boundary conditions given by where u is the x component and v is the y component of the Casson nanofluid velocity. U ∞ = bx is the free stream velocity of the Casson fluid, α is the thermal diffusivity of the fluid, μ nf is the effective dynamic viscosity of the Casson nanofluid, β is the non-Newtonian or Casson parameter, ρ nf is the effective density of the Casson nanofluid, k nf is the effective thermal conductivity of the nanofluid, k 1 is the permeability of the porous medium, F is the Forchheimer drag force coefficient, T is the temperature of the fluid, ðρC p Þ nf is the effective heat capacity of the Casson nanofluid, C p is the specific heat at a constant pressure of the fluid, μ f is the dynamic viscosity of the fluid, L is the slip length coefficient, k f is the thermal conductivity of the fluid, h f is the convective heat transfer coefficient, T f = T ∞ + nx 2 is the convective fluid temperature below the stretching/shrinking sheet, a is the constant of the linear stretching rate (s −1 ) of the Casson fluid, and b is the constant of the linear stretching/shrinking rate (s −1 ) of the wall (stretching for b > 0 and shrinking for b < 0). The parameters μ nf , ρ nf , k nf , and ðρC p Þ nf are defined by Tadesse et al. [37] as where ρ f , ρ s , k f , k s , ϕ and μ f are the density of the base fluid, density of the solid nanoparticle, base fluid thermal conductivity, nanoparticle thermal conductivity, nanoparticle volume fraction, and base fluid dynamic viscosity, respectively. Note that n = 3/Ψ, where Ψ, called the "sphericity," is defined as the ratio of the surface area of the sphere to that of the particle for the same volume. For spherical particles, Ψ = 1, and for the cylinders, Ψ = 0:5. This study considers the copper particle and is spherical in shape, so that n = 3, according to Hamilton and Crosser [40]. The thermophysical properties of nanoparticles and base fluid at T = 300K are given in Table 1, according to Tshivhi and Makinde [28] and Shaw et al. [34]. The above Equations (4)-(6) represent the governing equations of the nanofluid flow when β ⟶ ∞ and the Casson fluid flow when β ≠ ∞ and ϕ = 0. Equations (3)-(5) can be transformed to the dimensionless form by using the nondimensional variables below to obtain the similarity solutions by Tadesse et al. [37].
where the stream function ψ = x ffiffiffiffiffiffiffi aν f p f ðηÞ is written as a velocity component as

Computational and Mathematical Methods
Since the continuity equation is satisfied by (9) and (10), automatically, Equations (4) and (5) are converted into the following nondimensional form: With the dimensionless form of the boundary conditions where η is the similarity variable, λ is the velocity ratio (stretching/shrinking) parameter where λ > 0 for stretching and λ < 0 for shrinking of the sheet, Da is the Darcy number (porous media parameter), F r is the Forchheimer (second order porous resistance) parameter, Pr is the Prandtl number, Ec is the Eckert number, S is the constant mass flux parameter where S > 0 for suction and S < 0 for injection of the fluid, δ is the velocity slip parameter, and Bi is the Biot number (convective parameter). These dimensionless parameters and the variables A 1 , A 2 , A 3 , and A 4 quantities are defined as Note that in our discussion, F r , Da, Pr, Ec, and Bi measure the pressure drop caused by fluid-solid interactions to that of viscous and inertia resistance ratio, the relative effect of the permeability of the porous medium versus its cross-sectional area, momentum diffusivity to thermal diffusivity ratio, kinetic energy of the flow to heat dissipation potential (enthalpy difference) ratio across the thermal boundary layer, and internal thermal resistance at the surface of the sheet to the boundary layer thermal resistance ratio, respectively.
Moreover, according to [41], the pressure gradient in the flow due to porous medium ∇P is given by where q is the velocity vector, k 1 is the permeability of the porous medium m 2 , C F = 11/20ð1 − 11/20ðd/D e ÞÞ is a dimensionless form-drag constant, d is the diameter of spheres of the porous medium, and D e = 2wh/w + h is the equivalent diameter of the bed (defined in terms of the height h and width w of the bed). Thus, for our case, putting F ≡ ρC F (kgm -3 ), we obtain Fr which is dimensionless. The wall skin friction τ w and heat flux q w are computed as where Re x = xU ∞ /ν f is the local Reynold number. To compute the heat transfer enhancement (HTE) of the nanoparticles, we use the following formula: Fluid flow models of such kind could have dual solutions based on the physical parameters within the problem from the numerical result obtained. Therefore, by making a stability analysis, we determine the solution which is stable and physically practicable. We perform this analysis mathematically to validate the real solution among all the others. Thus, to employ the stability analysis, Equations (4) and (5) should be rewritten in unsteady (time dependent) case, according to Merkin [42]. Thus, we have where here t is time. (19) and (20) are transformed as follows:

Now, unsteady Equations
where τ is the nondimensional time variable. Using (21) in (19) and (20), we have with the boundary conditions In order to test the stability of solutions of f ðηÞ = f 0 ðηÞ and θðηÞ = θ 0 ðηÞ that satisfy boundary value problems (11)-(13), we have where ε is an unknown eigenvalue parameter (a small disturbance of growth or decay) that provides an infinite set of the 6 Computational and Mathematical Methods eigenvalues ε < ε 1 < ε 2 < ε 3 < …, and Fðη, τÞ and Gðη, τÞ are small relative to f 0 ðηÞ and θ 0 ðηÞ, respectively. The following linearized problem will be obtained by substituting (25) into (22)- (24).
A 3 Pr subjected to the boundary conditions Following Weidman et al. [43], the initial growth or decay of the solution (25) can be identified by obtaining the steady state solution setting τ = 0 so that F = F 0 ðηÞ and G = G 0 ðηÞ in Equations (26)-(28), where 0 < F 0 ðηÞ ≪ 1 and 0 < G 0 ðηÞ ≪ 1. The stability of the solution depends upon the sign of the smallest eigenvalue ε. If the value of ε 1 is positive, that shows the flow is stable, and there is an initial decay. Conversely, if the value of ε is negative, that shows the flow is unstable and illustrates an initial growth of disturbance. The linearized eigenvalue problem is given by subjecting to boundary conditions To obtain the smallest eigenvalues, there is a need to relax one of the boundary conditions in the form of the initial condition as suggested by Harris et al. [32]. In this problem, F 0 ′ ð∞Þ ⟶ 0 has been relaxed in the initial form as F 0 ″ ð0Þ = 1. The smallest negative eigenvalues point out the initial development of the disturbance, and the solution of the flow is unstable. On the other hand, the smallest positive related eigenvalue value shows the fluid flow is stable and physically realizable. Thus, the modified boundary conditions (31) become Equations (29) and (30) of linear eigenvalue problem are solved with the boundary conditions (31).

Numerical Method
To solve (11) and (12) along with the boundary condition (13) numerically, we use the method of fourth-fifth order using the shooting method with the MAPLE software. In order to apply the solver, first, we reduce the system as a set of equivalent ordinary differential equations of first order using the substitutions yð1Þ = f and yð4Þ = θ as below For the boundary conditions (13), To do stability analysis, we follow the same procedures. New substitutions are introduced to rewrite Equations (29) and (30) and the boundary conditions (32) into first order ordinary differential equations, by letting yð1Þ = F 0 , yð4Þ = G 0 , zð1Þ = f 0 , and zð4Þ = θ 0 .
For the boundary conditions, To determine the unknown initial conditions k 1 , k 2 , l 1 , and l 2 , we shoot them for arbitrary slope so that the solution of the system of ODEs satisfies the boundary conditions at ∞, and its accuracy is checked by comparing the calculated quantities with the provided end points. After obtaining these values, we apply the fourth-fifth order Runge-Kutta-Fehlberg technique to solve a system of first-order ODEs in (33) with boundary conditions (34) and determine ε from (35). To get the dual solutions, we take different initial approximates for the values of k 1 , k 2 ,where all profiles asymptotically satisfy the ∞ boundary conditions.

Results and Discussion
In this study, the results of the combined effects of velocity ratio (stretching/shrinking) parameter λ, Casson's parameter (factor) β, Darcy's number Da, Forchheimer's (porous medium inertia) parameter F r , suction/injection parameter In addition, the occurrence of two solutions for certain ranges of parameter variations is demonstrated for the coefficient of skin friction C f and Nusselt number (rate of heat transfer) Nu x in graphs and/or tables for diverse numerical quantities of parameters. The existence of dual solutions due to shrinking surface for certain ranges of parameter variations is shown for the skin friction coefficient using graphs and tables, and also, the variations of Nusselt's numbers are demonstrated in graphs for different values of change of parameters. The governing nonlinear ordinary differential equations, Equations (11) and (12), with the boundary conditions (13) are solved numerically using shooting techniques built-in Maple2018 software. The MAPLE solver has been used widely by many researchers to solve the boundary value problem (BVP), and this solver is coded with a finite difference in fourth order accuracy level. To validate this method, we compare the computational results obtained with that of the coefficient of skin friction of preceding works given by Nazar et al. [24], Jumana et al. [27], and Tadesse et al. [37], as presented in Table 2. It is observed that there is a nice agreement, and hence, the method is proper for tackling problems of the current study. The velocity and temperature profiles as well as skin friction and Nusselt's number are graphically presented. That means the interval of the shrinking parameter λ for which the similarity solution exists widens as suction and slipperiness of the surface parameters increase, whereas it diminishes as the Casson factor, nanoparticle volume fraction, porous medium, and porous medium inertia parameters increase. On the other side of this critical value λ c , no similar solutions exist because of the boundary layer separation from the surface of the sheet, and it is not possible to get the solution using the boundary layer approximations. In Figure 2(a), we see that for the upper branch solution for the shrinking parameter λ, the value of the skin friction coefficient (surface drag force) lessened with a rise in nanoparticle volume fraction ϕ and reversed as it approached the critical value of the shrinking parameter. In reality, the rising nanoparticle volume fraction means that the nanoparticles and the base fluid collide with each other, which raises the motion of the nanofluid; as a result, the thickness of the momentum boundary layer diminishes and increases the surface drag force. This result agrees with the computed and tabulated values in Table 3. For λ = 1, we always obtain f ″ ð0Þ = 0 for all values of ϕ since the fluid velocity is the same as velocity of the stretching/shrinking surface of the sheet. In Figure 2(b), we see that the skin friction coefficient diminishes as the Casson factor β increases for the upper branch solution. This is due to an increment in β, which means that the fluid loses its non-Newtonian behavior and acts like a Newtonian type, and hence, its velocity enhances, resulting in the reduction of shear stress.  Table 4. Figure 5 demonstrates the presence of a dual solution using the graph of local Nusselt's number Nu x , by taking particular numerical values of viscous dissipation parameter Ec and convective heating parameter Bi, to observe their impact on the process of heat transfer and the intervals of dual solution existence, for the shrinking sheet. It is noted from the model that the energy and momentum equations are coupled and hence the Nusselt number characterizes a dual solution for λ c < λ < 0 in the case of a shrinking sheet surface as displayed for particular governing parameters in figures. In Figure 5(a), it is observed that as viscous dissipation parameter Ec increases, the local Nusselt number enhances for the upper branch solution and drops for lower branch solutions, where the upper branch solution is the one which is stable and physically realizable, while the lower branch solution is not. Thus, minimum viscous dissipation (minimum Ec) is an acceptable sign for the rate of heat transfer enhancement of the nanofluid as the literature reveals. The interval of solution widens as the viscous dissipation Ec escalates. Again, Figure 5(b) reveals that as convective heating parameter Bi increases, the local Nusselt number rises for both the upper and lower     shrinking of the corresponding boundary layer thickness. This figure again describes how the porous medium inertia parameter Fr engulfs the velocity profile f ′ ðηÞ within the flow regime. Thus, the increment in porous media inertial resistance looms the velocity profile and the related reduced boundary layer thickness. Figure 7(a) depicts that intensifying the values of Cu-H 2 O volume fraction ϕ increases the velocity of the flow. That is, the higher ϕ augments the velocity profile due to the spherical-shape of Cu nanoparticle. Physically, the higher ϕ enhances the motion of the spherical-shaped the Cu nanoparticles which diminishes the boundary layer thickness and raises the velocity profile of the Cu-H 2 O-Casson nanofluid flow. The same figure also shows that improving the value of suction S augments the velocity of the flow, which is because suction is dominated by slipperiness of the surface and reduces the drag on the surface to control the separation of the boundary layer. As the slipperiness parameter δ rises, the velocity profile escalates, whereas its boundary layer thickness diminishes as demonstrated in Figure 7    14 Computational and Mathematical Methods due to the porous medium parameter Da is demonstrated in Figure 8(b). The temperature profile falls due to an increment in Da within a few boundary layers of fluid from the surface of the sheet (η < 1:13), and beyond (η > 1:13) the temperature profile switches trend due to an increment in Da. This figure also displays the impact of the nonlinear porous medium inertia parameter Fr on the fluid temperature profile. It is observed that an increment in the porous medium inertia parameter produces the falling fluid temperature profile and the thermal boundary layer. By looking at Figure 9(a), to observe how the temperature profile behaves in dealing with different Cu-H 2 O nanoparticle volume fractions ϕ, and it is seen that as ϕ increases, the temperature profile θðηÞ and the corresponding thermal boundary layer also increase. Since the Cu nanoparticle has high thermal conductivity, physically, it means increasing nanoparticle volume fraction ϕ upsurges the temperature profile and heat transfer rate. That is, due to enhanced thermal conductivity, Cu nanoparticles heighten the thermal enhancement of the Casson nanofluid past the permeable surface. The same figure reveals that augmented fluid suction parameter S produces both a diminutive temperature profile and its thermal boundary layer thickness. As a result, when impacts of viscous dissipation overrun, for augmenting suction parameter S, the thermal boundary layer thickness diminishes. The distribution of the temperature field gets more consistent within the boundary layer of flow because of the impact of suction. Figure 9(b) illustrates the impact of viscous dissipation (Eckert's number Ec) on temperature profile. There is an indication that an increase in the viscous dissipation parameter is highly beneficial to the growth of temperature profile in the flow region. We know that viscosity absorbs kinetic energy from the fluid motion and changes it into internal energy that improves heating the fluid flow, hence increasing the temperature profile. Moreover, the same figure demonstrates that augmenting the convective heating (increment in Biot number Bi) enhances the temperature profile of Cu-H 2 O-Casson nanofluid flow. That means, the coefficient of heat transfer caused by the hot fluid beneath the sheet is directly associated with the convective heating parameter Bi.

Existence of Dual Solutions due to Shrinking Sheet. In
That is, behavior of temperature profiles due to convective heating is higher near the boundary layer, implicating that the fluid temperature adjacent to the sheet surface surpasses that of the free stream fluid flow. Thus, the increasing convective heating parameter Biresulted in an increment in the convective temperature profile of the fluid flow. In Figure 10(a), it is observed that increasing sheet surface slipperiness δ reduces both the temperature profile and the thermal boundary layer thickness of the Cu-H 2 O-Casson nanofluid flow.

4.4.
Rate of Heat Transfer. Figure 10(b) reveals that the intensification in the values of the Casson factor β, porous media parameter Da, and porous media inertia parameter Fr resulted in the enhancement of the local Nusselt's number for the stretching sheet. An increment in the Casson factor β raises the fluid motion due to slipperiness and drops the thermal profile, resulting in enhanced heat transfer rate, whereas as the porous media parameter increases, both the fluid motion and the thermal profile diminish, which also enhances the Nusselt number (rate of heat transfer) rate at the stretching sheets surface. Figure 11(a) depicts that the rise in the values of nanoparticle volume fraction ϕ resulted in the rise of the Nusselt number Nu, which in reality means that the rate of heat transfer at the surface of the sheet enhances as nanoparticle volume fraction rises. In other words, the addition of the Cu-H 2 O nanoparticle volume fraction produces a rise in the temperature gradient at the sheet surface and results in the rate of heat transfer from the surface of the sheet to the fluid. Again, it can be concluded that the heat transfer properties of the Casson fluid becomes better because of the addition of the Cu-H 2 O nanoparticles into it. The same figure also revealed that rising in the values of stretching parameter λ resulted in the ascending of the Nusselt number Nu, which means as stretching of the surface increases, the heat transfer at the surface rises. Again, the stepping-up in the values of slip parameter δ reduces the Nusselt number Nu, and hence, the heat transfer rate drops, which is due to improvement in slipperiness of the sheet surface reducing adhesion of the fluid to it, as we see from the same figure. Figure 11(b) explains the impact of suction parameter S, Eckert's number Ec, and Biot's number Bi on the Nusselt number for stretching sheet surfaces. From this figure, it is observed that an increment in suction parameter S boosts the processes of heat transfer rate from the stretching sheet surface, which is due to an increase in the suction rate that raises the Cu-H 2 O nanoparticle accumulation on the stretching sheet surface that resulted in an enhancement in heat transfer rate (Nusselt number Nu) of the fluid. In other words, the minimum temperature gradient at the sheet surface is closely linked with an inverse relation to the heat transfer rate from the surface of the sheet. This shows that higher suction parameter S resulted in a considerably higher heat transfer rate. This figure also illustrates the impact of viscous dissipation on the management of heat transfer from the surface of the sheet. It is observed that the viscous dissipation parameter Ec raises the heat transfer rate falls, which is due to the slipperiness of the surface and the addition of nanoparticles, and moreover, the flow motion increases, and hence, the temperature gradient drops, which, in turn, results in an enhanced rate of heat transfer. Due to convective heating of the sheet surface, it is seen that the Nusselt number (heat transfer rate) rises as convective parameter Bi augments, as we see from the same figure. 4.5. Heat Transfer Enhancement. The chart in Figure 12 demonstrates 4.6. Numerical Analysis of Stability Test. From the computational results of this problem, the dual solution exists for some interval of λ. Temporal stability analysis is made to determine stable solutions among different solutions appear due to a shrinking sheet surface. As detailed in Table 3, for different values of Cu-H 2 O nanoparticle volume fraction ϕ and velocity ratio λ, the smallest eigenvalue ε is calculated for the temporary change of small perturbations regarding the basic steady flow, where the values of β = 10, Da = Fr = S = δ = 0:1 are fixed. Similar calculations were done for fixed ϕ = 0:1 as presented in Table 4. It is noted that since the value of the smallest eigenvalue ε is positive, the upper branch solution for shrinking sheet surface is hydrodynamically temporally stable and, hence, physically realizable. Clearly, for the lower branch solution, the negative value of ε revealed that it is unstable and physically unrealistic. In addition, ε > 0 demonstrates the rate of declination of small disturbance in the upper branch solution, whereas ε < 0 or the lower branch solution shows the increment of the disturbances.

Conclusions
A numerical investigation of boundary layer flow of Cu-H 2 O-Casson nanofluid past a slippery stretching/shrinking sheet through a Darcy-Forchheimer porous medium has worked out to demonstrate the overall impacts of Casson's factor, viscous dissipation, suction/injection, convective heating, and porous medium resistance parameters. Using similarity transformations, the modeled boundary layer equations are transformed into a system of ODEs, and the MAPLE software package is used for the computation of the numerical solutions. Stability analysis has been done to identify stable and physically reliable solutions subjected to small perturbations. The effects of various parameters on the dimensionless velocity and temperature profiles, coefficient of skin friction, the rate of heat transfer, and heat transfer enhancement are obtained numerically and presented in graphs, tables, charts. The following findings are derived from the discussions: (i) There is a critical value of shrinking parameter λ c below which no real, dual, or unique solutions exist. The critical value jλ c j widens with increase in suction, slipperiness of surface, and viscous dissipation parameters, whereas it diminishes with increase in Cu-H 2 O nanoparticle volume fraction, Casson's factor, porous medium, and porous medium inertia parameters (ii) Temporal stability analysis provided the smallest eigenvalue ε, which revealed that only the upper branch solution is stable and physically realizable, whereas the lower branch solution is unstable and not realistic (iii) The skin friction (wall shear stress) at the sheet surface increases as nanoparticle volume fraction and suction parameters escalate but reduces with increment in Casson's factor, slipperiness of the sheet, porous medium, and porous medium inertia parameters (iv) The velocity profile augments with the increment in nanoparticle volume fraction, velocity ratio, and suction and slipperiness parameters. However, it decreases with increasing values of the Casson factor, porous medium, and porous medium inertia parameters for the upper branch solution (v) Widening of the momentum boundary layer thickness is observed with enhancement in porous medium and porous medium inertia parameters, whereas it diminishes with increment in velocity ratio, suction, slipperiness of surface, Casson's factor, and nanoparticle volume fraction parameters (vi) The temperature profile and the thermal boundary layer thickness rise with incremental nanoparticle volume fraction, Casson's factor, viscous dissipation, and convective heating, whereas they decrease with increase in velocity ratio, suction, slipperiness, porous medium, and porous medium inertia parameters (vii) The rate of heat transfer intensifies with the increment in Casson's factor, porous medium, porous medium inertia parameters, suction, velocity ratio, nanoparticle volume fraction, and convective heating parameters but drops with increase in slipperiness of the surface and viscous dissipation parameters