Unsteady Three-Dimensional MHD Non-Axisymmetric Homann Stagnation Point Flow of a Hybrid Nanoﬂuid with Stability Analysis

: The hybrid nanoﬂuid under the inﬂuence of magnetohydrodynamics (MHD) is a new interest in the industrial sector due to its applications, such as in solar water heating and scraped surface heat exchangers. Thus, the present study accentuates the analysis of an unsteady three-dimensional MHD non-axisymmetric Homann stagnation point ﬂow of a hybrid Al 2 O 3 -Cu / H 2 O nanoﬂuid with stability analysis. By employing suitable similarity transformations, the governing mathematical model in the form of the partial di ﬀ erential equations are simpliﬁed into a system of ordinary di ﬀ erential equations. The simpliﬁed mathematical model is then solved numerically by the Matlab solver bvp4c function. This solving approach was proﬁcient in generating more than one solution when good initial guesses were provided. The numerical results presented signiﬁcant inﬂuences on the rate of heat transfer and ﬂuid ﬂow characteristics of a hybrid nanoﬂuid. The rate of heat transfer and the trend of the skin friction coe ﬃ cient improve with the increment of the nanoparticles’ concentration and the magnetic parameter; however, they deteriorate when the unsteadiness parameter increases. In contrast, the ratio of the escalation of the ambient ﬂuid strain rate to the plate was able to adjourn the boundary layer separation. The dual solutions (ﬁrst and second solutions) are obtainable when the surface of the sheet shrunk. A stability analysis is carried out to justify the stability of the dual solutions, and hence the ﬁrst solution is seen as physically reliable and stable, while the second solution is unstable.


Introduction
The stagnation point flow has attracted vast attention from many researchers because of its broad applications in both industrial and scientific applications.Some of the real-world applications of the stagnation point flow lie in the polymer industry, extrusion processes, plane counter-jets, and numerous forms of hydrodynamic modelling in engineering uses ( [1][2][3]).An exact solution of the steady two-dimensional stagnation-point flow towards a solid surface in moving fluid was first discovered by Hiemenz [4] in 1911.In his particular study, the Navier-Stokes equations are reduced to non-linear ordinary differential equations by using a similarity transformation.The remarkable work done by Hiemenz [4] was extended by Homann [5], who started the classical work of three-dimensional stagnation point flow for the axisymmetric case.Meanwhile, the flow in the neighbourhood of a particular stagnation point on a surface was explored by Howarth [6], focusing on the non-axisymmetric three-dimensional flow near the stagnation region.Fast forward to 1961, the work of Howarth [6] was criticized by Davey [7], who indicated a mistake in Howarth's paper exposing that the results in the region −1 ≤ c ≤ 0 are unable to be achieved from those discovered for 0 ≤ c ≤ 1 as reported in the study.In conjunction with these findings, Davey and Schofield [8] initiated the study of these saddle point solutions and justified the existence of the non-uniqueness solution.In another study, Weidman [9] modified the Homann's axisymmetric outer potential stagnation-point flow for non-axisymmetric stagnation flow of the strain rate.The study revealed a new clan of asymmetric viscous stagnation point flows liable on the shear rate ratio, γ = b/a where −∞ ≤ γ ≤ ∞, a is the strain rate and b is the shear rate.An analysis of unsteady heat transmission in non-axisymmetric Homann stagnation-point flows of a viscous fluid over a rigid plate was investigated by Mahapatra and Sidui [10], and recently, an investigation on the non-axisymmetric Homann stagnation-point flows of a viscoelastic fluid towards a fixed plate was conducted by Mahapatra and Sidui [11].
Ever since the evolution study of the stagnation point flow with the presence of dual solutions by Davey [7], various works concerning the stagnation point flow towards a shrinking sheet were introduced.Wang [12] considered two-dimensional stagnation point flow on a two-dimensional shrinking sheet and axisymmetric stagnation point flow on an axisymmetric shrinking sheet, while Mahapatra and Sidui [13] assessed unsteady heat transfer in non-axisymmetric Homann stagnation-point flow towards a stretching/shrinking sheet with stability analysis.The continuous effort was carried out by Khashi'ie et al. [14] who examined the three-dimensional non-axisymmetric Homann stagnation point flow and heat transfer past a stretching/shrinking sheet using hybrid nanofluid.Meanwhile, Zaimi and Ishak [15] scrutinized the slip effects on the stagnation point flow towards a stretching vertical sheet.Nevertheless, explorations on the stagnation point flow keep evolving in various ways and have been working still because of its importance in massive engineering applications and also in the magnetohydrodynamics (MHD) flow field.A comprehensive study of the literature on the related works was reviewed by [16][17][18][19].
A fluid that is heated by electric energy in the occurrence of a vigorous magnetic field, such as crystal growth in melting, is essential in the industrial sector.The interaction of electrical currents and magnetic fields generates the divergence of Lorentz forces during the movement of fluid.In accordance with this phenomenon, MHD describes the hydrodynamics of a conducting fluid in the presence of a magnetic field.The examinations of MHD flow are very significant due to its massive number of uses implicating the magnetic effect in industrial and engineering areas, such as MHD electricity generators, sterilization tools, magnetic resonance graphs, MHD flow meters, and also in granular insulation (see [20,21]).The goods of the end product depend immensely on the rate of cooling involved in these processes, managed by the application of the magnetic field and electrically conducting fluids.The study of MHD flow in the Newtonian fluid was first carried out by Pavlov [22], who investigated the magnetohydrodynamic flow of an impressible viscous fluid caused by deformation of a surface.Chakrabarti and Gupta [23] broadened the study of hydromagnetic flow and heat transfer over a stretching sheet, followed by Vajravelu [24] who widened the hydromagnetic flow study over a continuous, moving, porous flat surface.Andersson, in 1995, introduced an exact solution of the Navier-Stokes equations for magnetohydrodynamic flow [25], and Lok et al. [26] analyzed the MHD stagnation-point flow towards a shrinking sheet using the Keller-box method and proved the existence of multiple (dual) solutions for small values of the magnetic field parameter for the shrinking case.Recently, Almutairi et al. [27] studied the influence of second-order velocity slip on the MHD flow of a nanofluid in a porous medium by considering the homogeneous-heterogeneous reactions.On the other hand, the impact of nonlinear and temperature jump on non-Newtonian MHD nanofluid flow and heat transfer past a stretched thin sheet was examined by Zhu et al. [28] Over the last few decades, the researcher has experienced tremendous scholarly devotion to the study of heat transfer fluid.Recent demand for a high-efficiency refrigeration system and the ineffectiveness of traditional thermal conduction fluids encouraged analysts to discover another heat transfer fluid.Choi and Eastman [29] launched the exploration of nanofluids and illustrated the presence of suspended nanoparticles in a carrier fluid.This pioneering study led to the verdict of the colloidal suspension of intensely small-sized particles, for instance, carbon nanotubes, metals, oxides, and carbides, into the based fluid, which may ensure access to an advanced course of nanotechnology-based heat transfer media (see [30,31]).The eccentric features of nanofluids have gained great acknowledgement in various engineering, medical, and industrial applications like engine cooling, diesel generator efficiency, micro-manufacturing, solar water heating, cancer treatment, nuclear reactors, and diverse types of heat exchangers ( [32][33][34]).Due to the massive potential for the applications of nanofluids, Choi and Eastmen [29] developed a mathematical model of nanofluids, which allowed Buongiorno [31] to contribute to heat transfer analysis in nanofluids by introducing the non-homogeneous model for transport and heat transfer phenomena in nanofluids with turbulence applications.
Recently, an expansion of new engineered nanofluids was achieved by dispersing composite nanopowder or dissimilar nanoparticles with sizes between 1 and 10nm in the base fluid [29]; it is known as a hybrid nanofluid.The hybrid nanofluid is a modern technology fluid that may offer better heat transfer performance and thermal physical properties.The progress related to the preparation methods of hybrid nanofluids, thermo-physical properties of hybrid nanofluids, and current applications of hybrid nanofluids was published by Sarkar et al. [35] and Sidik et al. [36].In another study, Huminic and Huminic [37] highlighted the essential applications of hybrid nanofluids, such as in heat pipes, mini-channel heat sinks, plate heat exchangers, air conditioning systems, tubular heat exchangers, shell and tube heat exchangers, tube in a tube heat exchangers, and coiled heat exchangers.Turcu et al. [38], for the first time testified the hybrid nanocomposite particle synthesis, which consisted of two different hybrids, polypyrrole-carbon nanotube (PPY-CNT) nanocomposite and multi-walled carbon nanotube (MWCNT) on magnetic Fe 3 O 4 nanoparticles.In the following year, Yen et al. [39] inspected the effect of hybrid nanofluids in channel flow numerically.Devi and Devi [40] analysed the problem of hydromagnetic hybrid nanofluid (Cu-Al 2 O 3 /water) flow on a permeable stretching sheet subject to Newtonian heating, and they continued the investigation to improve the heat transfer in hybrid nanofluid flow past a stretching sheet [41].Subsequently, Yousefi et al. [42] reviewed on the stagnation point flow of an aqueous titania-copper hybrid nanofluid toward a wavy cylinder.At the same time, Khashi'ie et al. [43] performed a numerical study on the heat transfer and boundary layer flow of axisymmetric hybrid nanofluids driven by a stretching/shrinking disc.A detailed documentary on the numerical study of hybrid nanofluid flow and heat transfer is reviewed by [44][45][46][47][48].
Many practical situations, such as a sudden stretching of the plate or temperature change of the plate, involved unsteady conditions of the heat transfer flow.Cai et al. [49] explained that the flow in the viscous boundary layer near the plate would slowly be enlarged if the surface was extended unexpectedly, and hence converted into a steady flow after a certain interval.Technically, we believe that the consideration of physical quantities related to time is crucial in mathematical modeling and analysis, which is acknowledged in the formulation of this research problem.
Motivated by the work by Mahapatra and Sidui [13], this study aims to inspect the unsteady MHD non-axisymmetric Homann stagnation point of a hybrid nanofluid in three-dimensional flow.The proposed hybrid nanofluid model is adapted from Devi and Devi [40] and Hayat and Nadeem [50], recognized by suspending varied nanoparticles, namely alumina (Al 2 O 3 ) and copper (Cu), in the base fluid (water).To the best of the authors' knowledge, no attempt has been made to examine the heat transfer and fluid flow of the hybrid nanofluid (Al 2 O 3 -Cu/H 2 O) considering the unsteady parameter in non-axisymmetric Homann stagnation point flow.This possibly will benefit future works on choosing a significant parameter to enhance the heat transfer performance in the modern industry.The novelty of this study can also be seen in the discovery of dual solutions and the execution of stability analysis.Ultimately, this research is highly claimed to be authentic and original.

Mathematical Model
Consider the unsteady three-dimensional MHD non-axisymmetric Homann stagnation point flow of a hybrid Al 2 O 3 -Cu/water nanofluid with a stretching/shrinking sheet on the x, y− plane where x, y and z are Cartesian coordinates with the z− axis measured in the horizontal direction and the axes x and y are in the plane z = 0 as illustrated in Figure 1, respectively.We assume that the constant surface temperature T w is stretched and shrunk in the x and y directions by the velocities u w = εcx 1+αt and v w = εcy 1+αt The uniform temperature is given by T ∞ and B 0 is introduced to the stretching/shrinking sheet in an orthogonal direction as a transverse uniform magnetic field.Meanwhile, the modified non-asymmetrically free streamflow along the x, y and z axes is described by ( [9]): (1) ( ) ( ) ) The velocity component in x − direction is given by u , while v is in y − direction.Next, the boundary conditions are: Here, a is the strain rate and b is the shear rate of stagnation point flow, correspondingly.By adapting the Tiwari and Das [30] nanofluid model, the continuity, momentum, and the energy equations of the hybrid nanofluid can be written as follows: The velocity component in x− direction is given by u, while v is in y− direction.Next, the boundary conditions are: Note that T is the hybrid nanofluid temperature, µ hn f is the dynamic viscosity of the hybrid nanofluid, ρ hn f is the density of the hybrid nanofluid, k hn f is the thermal conductivity of the hybrid nanofluid, (ρC p ) hn f is the heat capacity of the hybrid nanofluid, σ hn f is the electrical conductivity of the hybrid nanofluid, and the time-dependent of a transverse magnetic field is given by The hybrid nanofluids thermophysical properties are specified in Table 1, as demonstrated by Devi and Devi [41,51].At this point, φ is the volume fraction of nanoparticles, ρ f and ρ s are the base fluid density and hybrid nanoparticles, C p is the heat capacity, (ρC p ) f and (ρC p ) s represent capacitance heating of the base fluid and hybrid nanoparticles, and finally k f and k s are the thermal conductivities of the base fluid and hybrid nanoparticles, respectively.Meanwhile, the thermophysical properties of the fluid and nanoparticles for aluminum oxide, copper, and the base fluid (water) are given in Table 2.

Properties
Hybrid Nanofluid Heat capacity Dynamic viscosity Thermal conductivity , where, Now, pursuing Mahapatra and Sidui [10], the resulting similarity transformation is proposed to achieve the similarity solutions: where the prime denotes differentiation with respect to η.By substituting (7) into the steady-state Equations ( 2)-( 5), the following ordinary differential equations are obtained: 1 Pr with the boundary conditions (6) which are converted to: In the equations mentioned above, A = α/c represents the unsteadiness parameter, λ = a/c is a ratio of the surrounding fluid strain rate to the surface strain rate, γ = b/c is the surrounding fluid shear rate ratio to the strain rate of the sheet, M = σ hn f /σ f ρ hn f /ρ f B 2 0 c denotes the magnetic parameter and Pr = ν f /α f indicates the Prandtl number.The parameter of stretching/shrinking is meant by ε where ε > 0 determines the stretching sheet, while ε < 0 reflects the shrinking sheet.The related quantities of interest in this study are the skin friction coefficient, C f x and C f y along the x− and y− directions and the local Nusselt number Nu x , which is specified as where τ wx , τ wy are the shear stresses along the x−, y− axes and q w represents the heat flux, correspondingly.Such terms can be defined by By prompting Equations ( 7), ( 12) and ( 13), we get: where Re x =

Stability Analysis
The system of Equations ( 8)- (10) along with the boundary conditions (11), is capable of generating more than one solution and ultimately permits the requirement analysis of the flow to identify the reliable and feasible solution.Going through the outstanding work done by Merkin [53] and Merrill et al. [54] in the stability analysis, the application of an unstable form of the boundary layer problem is analyzed by using the time variable and a dimensionless time variable denoted by τ [55].Next, we consider the following new similarity variables: Using the similarity variables of Equation ( 15) into Equations ( 8)-( 10), the altered differential equations and the boundary conditions are as follows: 1 Pr Weidman et al. [55] highlighted that the stability of solutions is introduced by determining the system's decay or initial growth.This can be achieved via considering the following perturbation expressions of the primary flow f = f 0 (η), g = g 0 (η) and θ = θ 0 (η) with the resulting equation: where ω is an unknown parameter of the eigenvalue, F(η), G(η) and H(η) are comparatively slight to f 0 (η), g 0 (η) and θ 0 (η).Substituting (20) into Equations ( 16)-( 18) we attained the following system of equations: 1 Pr subject to the boundary conditions: The stability of the steady-state flow and heat transfer solutions f 0 (η), g 0 (η) and θ 0 (η) was implemented by setting τ = 0 with F = f 0 (η), G = g 0 (η) and H = θ 0 (η) in ( 21)- (24).Finally, the initial development or the solution decay of the solution (20) is identified.The value of ω is obtained by solving the following eigenvalue problem: 1 Pr along with the conditions: The range of possible eigenvalues can be determined by relaxing a boundary condition on F 0 (η) according to the previous research by Harris et al. [56].In this study, the condition F 0 (∞) → 0 is relaxed, and the linear eigenvalue problem ( 25)-( 28) is solved together with the new boundary condition F 0 (0) = 1 for a fixed value of ω 1 .The flow is considered unstable if the smallest eigenvalue ω 1 .is negative, which indicates an initial growth of disturbances occurred, while a positive value of the smallest eigenvalue signifies that the flow is physically achievable and stable.

Results and Discussion
The bvp4c function in Matlab is adopted to produce the results of the nonlinear system of ordinary differential Equations ( 8)-( 10) together with the boundary conditions (11).The relative error tolerance is set as 10 −10 to gain the results of the numerical outcomes and stability analysis.Table 3 presents that the average central processing unit (CPU) time required for computing each result in Table 4 is approximately 1.4 s, and the dual solutions were obtained by indicating different initial guesses which must satisfy the far-field boundary conditions (11) asymptotically.The numerical results are validated with the numerical results produced by Mahapatra and Sidui [10] and Nawaz and Hayat [57] by collating the values of the shear stress f (0) of an axisymmetric (γ = 0) stagnation point flow of a viscous fluid with the exclusion of magnetic field and the unsteady parameter, which is depicted in Table 4.It demonstrates excellent agreement with the previous literature; hence, the practicality and effectiveness of the bvp4c method are verified.The estimated relative error, ε r is also measured, and it shows that the calculated values of ε r are relatively small between the present and previous results.Figure 2 exhibits the variation of the wall shear stress parameter f (0) and g (0) towards the difference value of γ when φ 1 = φ 2 = M = A = 0, λ = 0.1, ε = 1.0 and Pr = 1.0 where the dotted lines correspond to the asymptotic behaviour of f (0) and g (0) which is in excellent agreement with Mahapatra and Sidui [10] who pursued a standard fourth-order Runge-Kutta integration technique in their study.This justifies the role of the bvp4c numerical technique as a dependable practice, and the present results are valid and correct.The dimensionless velocity profiles f (η) and g (η) for different values of λ are illustrated in Figures 3 and 4. The figures prove that the profiles of the velocity comply with the far-field boundary conditions of Equation ( 11) asymptotically.The maximum value of the velocity gradient with the lowest thickness of the momentum boundary layer is preserved for the largest value of λ.Notably, the distance of two adjacent profiles increases remarkably with the increased amount of λ in both the first and second solutions.Figures 5 and 6 portray the variations of f (0) and g (0) towards ε for different values of λ in hybrid nanofluids with the existence of the magnetic field and the influences of the unsteadiness parameter.An enrichment in the amounts of λ embarking on the augmentation of both λ + γ and λ − γ provide a significant effect to the surface shear stresses.The variations of f (0) are expected to be higher by the increasing value of λ + γ yet decrease when λ + γ is decreased, and the same trend is expected in g (0) with λ − γ.On the other hand, the effect of the nanoparticles volume fraction is observed in Figures 7-9, respectively.Surprisingly, the critical values of the various usage of the nanoparticles' volume fraction give no significant effect to the trend of the nanofluid (φ 1 = 0.02, φ 2 = 0.00) and hybrid nanofluid flow (φ 1 = 0.02, φ 2 = 0.002, 0.04).It is observed that the reduced skin friction coefficient in x and y directions, as presented in Figures 7 and 8, and the reduced local Nusselt number in Figure 9 upsurges with the rising in the nanoparticles' volume fraction.This is due to the fact that more kinetic energy is produced with a higher concentration of nanoparticles and thus, enhances the heat transfer of the fluid particles.This finding also corresponds with several present particle-laden direct numerical studies (DNS), which reveal that the roughness components appear to redistribute the energy and, therefore, decrease the overall large-scale near-wall anisotropy of the flow pattern (see Yuan and Piomelli [58] and Ghodke and Apte [59]).Table 4. Results of f (0) for specific values of λ when φ 1 = 0.0, φ 2 = 0.0, γ= 0, ε= 1, M = 0, A = 0, and Pr = 6.2.

λ
Present Result Mahapatra and Sidui [10] ε r =|     energy is produced with a higher concentration of nanoparticles and thus, enhances the heat transfer of the fluid particles.This finding also corresponds with several present particle-laden direct numerical studies (DNS), which reveal that the roughness components appear to redistribute the energy and, therefore, decrease the overall large-scale near-wall anisotropy of the flow pattern (see Yuan and Piomelli [58] and Ghodke and Apte [59]).Figures 10-12 exhibit the influence of the magnetic parameter on f (0), g (0) and −θ (0) which shows a prominent effect on the fluid flow of the shrinking sheet.The reduced skin friction coefficient in both the x− and y− directions rise and eventually increase the value of −θ (0) with the escalation of M due to the occurrence of the Lorentz force, which acts to retard the fluid flow.The Lorentz force creates resistance to the motion of the fluid particles and then consequently reduces the fluid velocity.The synchronism of the magnetic and electric field that occurred from the formation of the Lorentz force tends to slow down the fluid movement.The boundary layer becomes thinner in both directions as proven in Figures 13 and 14 as M increases due to the delayed flow, hence contributing to the increment of f (0) and g (0).On the other hand, Figure 15 advertises the variations of the non-dimensional temperature profiles for different values of M which tend to decrease in the first solution but show a reverse trend in the second solution.The figures clearly reveal that the thicknesses of the thermal boundary layers decrease with the increase of M Practically, by restricting the magnetic field intensity, the progression of the thermal boundary layers' thicknesses can be managed and, thus, be able to reduce the temperature profile distributions.The effect of the unsteadiness parameter A on f (0), g (0) and −θ (0) is depicted in Figures 16-18, respectively.Increasing the values of A results in a reduction of the skin friction coefficients in both directions, as illustrated in Figures 16 and 17, which consequently decreases the reduced local Nusselt number.The fact that the unsteadiness parameter affects the velocity and temperature profile is proven in Figures 19-21.The velocity of the fluid is in decline along the surface of the sheet due to the increasing value of the shear stress and subsequently shrinks the thickness of the momentum boundary layer nearby the wall, as displayed in Figures 19 and 20, accordingly.The increasing value of A decreases in the temperature profile of the hybrid nanofluid, as shown in Figure 21, which is understandable because the spaces between the molecules become higher in unsteady flow and proportionately decrease the temperature profile and improve the cooling rates of the fluid.Therefore, the unsteadiness parameter should be highlighted and well considered for practical purposes.Figure 21, which is understandable because the spaces between the molecules become higher in unsteady flow and proportionately decrease the temperature profile and improve the cooling rates of the fluid.Therefore, the unsteadiness parameter should be highlighted and well considered for practical purposes.Since the dual solutions noticeably exist, as illustrated in Figures 3-21, a stability analysis is performed to discover a significant solution between the first and second solutions.This process is achievable by clarifying the eigenvalue problems in Equations ( 25)-( 27) using bvp4c in Matlab, which produce an infinite set of ω 1 < ω 2 < ω 3 . . . . ... The stability of the flow is reliant on the smallest positive eigenvalue ω 1 since there exists an initial decay of disturbances.In contrast, there is an initial growth of perturbations if the smallest eigenvalue ω 1 is negative, which signifies that the solution is unstable.As depicted in Table 5, the lowest eigenvalues for the first solution are positive, while the second solution is negative.In conclusion, the first solution is stable and physically reliable, while the second solution is unstable and unreal.

Conclusions
The current analysis is devoted to examining the unsteady three-dimensional non-axisymmetric Homann stagnation point flow of alumina (Al 2 O 3 ) and copper (Cu) hybrid nanofluids in the presence of MHD.The governing partial differential equations are transformed into a system of ordinary differential equations by using a similarity transformation and appropriately solved by a bvp4c function in Matlab.An increment in λ may upsurge the velocity gradient and thus decline the momentum boundary layer thickness.A high concentration of the nanoparticle volume fraction speeds up the molecules' kinetic energy and then enhances the heat transfer process of the fluid particles.The increment in the intensity of the magnetic parameter M increases the local Nusselt number and the skin friction coefficient.

Figure 1 .Figure 1 .
Figure 1.Flow model and coordinate systems of the physical model.Note that T is the hybrid nanofluid temperature, hnf μ is the dynamic viscosity of the hybrid

ε
is the estimate percentage relative error between the present result, r and previous result, s .

Figure 4 .
Figure 4. Velocity profiles of g (η) for different values of λ.

Figure 10 .
Figure 10.Variations of f (0) for different values of M.

Figure 12
Figure 12.Variations of

Figure 12
Figure 12.Variations of

Figure 14 .
Figure 14.Velocity profiles of g (η) for different values of M.

Figure 16 .
Figure 16.Variations of f (0) for different values of A.

Figure 16 .
Figure 16.Variations of (0) f ′′ for different values of A .

Figure 18
Figure 18.Variations of

Figure 18
Figure 18.Variations of

Figure 20 .
Figure 20.Velocity profiles of g (η) for different values of A.

Table 3 .
Computational time to generate the values of

Table 4 .
Results of