MHD Mixed Convection of Hybrid Ferrofluid Flow over an Exponentially Stretching/Shrinking Surface with Heat Source/Sink and Velocity Slip

Abstract

This paper investigates the mixed convection flow related to ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4/\mathrm{water}$ hybrid ferrofluids on stagnation point over an exponentially stretching/shrinking surface with the influence of magnetohydrodynamic (MHD), velocity slip, and heat source/sink. The proposed system of differential equations is reduced using the similarity transformation procedure that is numerically solvable using MATLAB. Dual solutions are obtained for various governing parameters, in which the first solution is found to be in a stable state via the stability analysis. The CoFe₂O₄ nanoparticles volume fraction increases the heat transfer rate and the skin friction coefficient while delaying the separation of the boundary layer at the bifurcation point. Adding CoFe₂O₄ nanoparticles in hybrid ferrofluids gives a better heat transfer rate than that obtained with ferrofluids. The presence of a magnetic field enhances the fluid flow velocity. The increased strengths of the heat sink and stretching parameters give better results on the heat transfer, while the results are reversed for the heat source and shrinking parameters. The presence of velocity slip does influence the skin friction and the fluid flow.

1. Introduction

Due to its industrial significance, nanotechnology has become more prevalent in various industries. Nanotechnology involves the use of nanoparticles and a base fluid to transfer heat, which has become a new category of heat transfer fluids. Choi and Eastman [1] introduced the term nanofluids as a way to introduce an entirely new type of fluid. Saidur et al. [2] also demonstrated that nanofluids have greater thermophysical properties and heat transfer efficiency than traditional fluids (water, glycerol). Moreover, Sidik et al. [3] and Sidik et al. [4] reviewed the application of nanofluids. For those reasons, numerous investigations have been accomplished to study boundary layer flow and heat transfer in nanofluids [5,6,7,8]. In contrast, magnetic nanofluid or ferrofluid is a suspension of magnetic nanoparticles held together by a colloidal relationship (magnetite, hematite, or cobalt ferrite) in a non-magnetic base fluid (ethylene glycol, oil, or water) [9]. In heat transfer problems, research into ferrofluids has received much attention because of the multitude of applications for ferrofluids in nanotechnology industries, such as heat exchangers, coolants in thermal management devices, and operations related to the transfer of heat, including boiling and improvised cooling functions of loudspeakers [10,11]. More significantly, the use of ferrofluids in biomedical fields has been highlighted since they make it easier to diagnose and treat diseases including rheumatoid arthritis, cancer, and angiocardiopathy [12]. Ferrofluid research is, therefore, thought to have a bright future and hold the potential to further transdisciplinary studies in biology, chemistry, physics, and material science. Due to this reason, several academicians and researchers have conducted experimental investigations on the behavior of ferrofluids toward heat transfer problems [13,14,15,16,17].

Recently, as a new ferrofluid category, hybrid ferrofluids have been investigated in numerous studies by researchers to assess their relationships with heat transfer and thermal efficiency. Hybrid ferrofluids are substances composed of two or more magnetic nanoparticles in a base fluid. Hybrid ferrofluids have been studied in numerous heat exchange processes and are used for a wide range of purposes. Many researchers are interested in using hybrid ferrofluids to address real-world heat transfer problems, with some of these noteworthy applications including damping, dynamic sealing, heat dissipation, etc. These applications have attracted the attention of many researchers regarding the use of hybrid ferrofluids. However, an examination of the published literature indicates that there are still very few papers about hybrid ferrofluids’ flow and heat transfer. The liquid film flow and heat transfer in a magnetohydrodynamic (MHD) hybrid ferrofluid in the presence of radiation and an irregular heat source/sink were studied by Kumar et al. [18]. They noticed that compared to ferrofluid, a hybrid ferrofluid delivers a better rate of heat transmission. Following these, a continued study by Tlili et al. [19] was conducted on an MHD hybrid ferrofluid with the influence of asymmetrical heat rise/fall. Furthermore, Manh et al. [20] examined the radiation effect in the heat transfer of an MHD hybrid ferrofluid immersed in porous media. An experimental procedure with aqueous hybrid ferrofluid was performed by Giwa et al. [21]. The authors considered the effect of uniform magnetic induction in a rectangular cavity and obtained positive results for the heat performance of hybrid ferrofluid compared to mono-nanoparticle. Recent work by Anuar et al. [22] studied the stagnation flow of a hybrid ferrofluid past an exponentially stretching/shrinking sheet under the influence of external heat and a magnetic field. They discovered from their research that an increase in the volume percentage of ${\mathrm{CoFe}}_2{\mathrm{O}}_4$ nanoparticles tended to postpone the boundary layer separation.

Studies regarding the flow characteristics and heat transfer of mixed convection have been carried out due to the great practical interest. Mixed convection is a combination of natural (free) and forced convection. It occurs when free and forced convection mechanisms are at work at a given time to contribute to flow and heat transfer. For example, Sparrow et al. [23] applied the aiding and opposing flow in mixed convection to study boundary layer flow. Afterwards, researchers’ interest in mixed convection rose, providing some notable mention [24,25,26,27]. Since then, with the hope of greatly enhancing aspects of boundary layer and heat transfer, researchers have already begun considering MHD in mixed convective flow. Moreover, the MHD term is essential in checking the heat transfer in creating value products (MHD generator, nuclear reactor, MHD power generation), as it relies upon the heat observing variable [28]. By accounting for the effect of the magnetic field in the mixed convection flow, Khan and Rasheed [29] looked into the heat transfer performance of nanofluid in MHD mixed convection flow with thermal radiation. They discovered that the presence of MHD reduced the skin friction and slowed the fluid flow velocity due to the resistance from the Lorentz force in flow. Furthermore, Wahid et al. [30] conducted a study on the effects of thermal radiation on mixed convection flow of hybrid nanofluid in the presence of MHD. They noticed from their study that as thermal radiation and magnetic values rose, the heat transfer performance was enhanced. Over the years, the study of boundary layer and heat transfer in MHD mixed convection flow has been an interesting subject explored among researchers and provides some quality mentions [31,32,33,34].

The reliability of the heat control factor-like heat source/sink is important to ensure end-product quality in the manufacturing industry, such as for nuclear reactors, rocket thrust systems, electronic chips, semi-conductor wafers, and many more. For example, Sharma et al. [35] explored the impact of heat source/sink on stagnation-point flow with the surface’s condition in vertical stretching. Then, Jamaludin et al. [36] further investigated the influence of heat source/sink over a permeable stretching/shrinking surface with thermal radiation and magnetic field. Meanwhile, Jamaludin et al. [37] also addressed the issue of mixed convection stagnation-point flow for hybrid nanofluid, $\mathrm{Cu}-{\mathrm{Al}}_2{\mathrm{O}}_3/\mathrm{water}$ in the presence of a heat source/sink and magnetic field over a permeable stretching/shrinking surface. Recently, Armaghani et al. [38] applied the L-shape cavity for testing the influence of a heat source/sink under the condition of mixed convection flow in hybrid nanofluid ${\mathrm{Al}}_2{\mathrm{O}}_3-\mathrm{Cu}/\mathrm{water}$, while Zainal et al. [39] considered the unsteady state of rear stagnation-point flow in a hybrid nanofluid with the influence of heat source/sink and magnetic field. Notably, there are some interesting articles that consider the heat source/sink conditions [40,41,42].

In the manufacturing industry (hot rolling, extrusion, and wire drawing), heat transfer over a moving surface to a quiescent ambient medium has already drawn much interest from most manufacturing sectors, since the temperature is higher and energy transfer tends to occur at the moving surface compared to ambient medium [43]. For example, Sakiadis [44] was the first to study the condition of moving surfaces on a boundary layer by considering a uniform speed over a continuous surface and solving it using an integral approach. A continuing study on [44] was conducted by Crane [45] using a two-dimensional boundary layer form on a linear stretching sheet. The first study of heat and mass transfer across an exponentially stretching surface with exponential velocity was performed by Magyari and Keller [46]. Then, the work of [46] was extended by Elbashbeshy [47] with the presence of a suction effect. For a new class of nanofluids, hybrid nanofluids, Anuar et al. [48] explored the existence of a dual solution to stagnation-point flow over a sheet that was exponentially stretching/shrinking with the condition of suction/injection. Of note, there are some interesting articles that consider the moving surface in conditions of exponential stretching/shrinking [49,50,51,52].

The appearance of slip has already been recognized as important over the past decades in the manufacturing industry, since the condition of no-slip boundary in fluid mechanics is not always applicable for several cases in complex fluids. Undoubtedly, the slip effect tends to appear in the system, as mentioned by Hatzikiriakos [53]. Complex fluids such as emulsions, suspensions, foams, and polymer solutions tend to form a wall slip, which exhibits a boundary slip at the wall. This phenomenon has importance in technological applications, for instance, the polishing of artificial heart valves and internal cavities [54,55]. The first studies on linear slip boundary conditions were conducted by Navier [56] and Maxwell [57]. Then, a critical analysis with a comprehensive discussion and argument between no-slip and slip boundary conditions was performed by Rao and Rajagopal [58] in their report. Later, Anuar et al. [59] investigated the heat transfer of carbon nanotubes with the influence of velocity slip and suction effect on mixed convection flow past an exponentially stretching/shrinking surface. It is worth noting that there are some interesting articles on the influence of velocity slips [60,61,62].

Finding dual or non-unique solutions is crucial for studying boundary layer theoretical models, because these solutions can highlight disrupted physical processes past moving surfaces. Particularly in the regulated arrangement, the interrupted physical process expresses the non-physically realizable condition. The disturbed physical process may also briefly appear in practical uses. In order to maintain the productivity of the particular production process, it is crucial to locate dual solutions. In order to make it easier to acquire the dual solutions, the boundary layer model must be solved numerically using the best approach. Convection fluid flow and heat transfer problems can be solved using a variety of methods, such as the Poisson–Boltzmann (PB) theory [63] and the traditional density functional theory (DFT) method. The shooting method was known to produce multiple numerical solutions in the past. The shooting technique, however, loses effectiveness when it comes to making early predictions. Those first guesses must be accurate and must be relatively near to the desired outcome or else dual solutions cannot be found in any other case. This is undoubtedly a time-consuming affair. As Ha [64] mentioned, when tackling a problem sensitive to the initial conditions, the shooting method’s other limitation becomes clear. The bvp4c function in MATLAB was used to solve the current problem, due to its effectiveness in generating several numerical solutions, and regardless of how complex the mathematical model is constructed, finding the non-uniqueness solution is possible. The adaptable bvp4c approach works well for solving a mathematical model with unknowable parameters. The reduced boundary layer mathematical models are solved using the bvp4c function, which is superior to other built-in numerical methods since it has a continuation code [65] that speeds up processing.

From the author’s knowledge and understanding of the recent literature, no previous research includes the mixed convection towards stagnation-point flow in the boundary layer and heat transfer problem in hybrid ferrofluids past an exponentially stretching/shrinking surface. As inspired by the literature mentioned above, the mathematical models will include the mixed convection towards stagnation-point flow in an ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4/\mathrm{water}$ hybrid ferrofluid for this present work. Other effects, for instance, heat source/sink, magnetic field, and velocity slip past an exponentially stretching/shrinking surface, will be considered in this problem. The mathematical models based on Tiwari and Das’s model [25] will be analyzed. Consequently, the stability analysis will be performed on dual solutions initiated by Merkin [66] and Weidman and Turner [67]. The method of obtaining stability between dual solutions plays a crucial role. It provides a procedure to determine which solution is stable, as already shown by other researchers [22,37,50] in their excellent papers on stability analysis.

2. Mathematical Formulation

A steady, incompressible, and two-dimensional magnetohydrodynamic (MHD) mixed convection flow towards stagnation point over a surface that is stretching/shrinking exponentially in a water-based hybrid ferrofluid (${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4/\mathrm{water}$) is considered in this present work, as displayed in Figure 1. The $x$- and $y$-axes in this case are the Cartesian axes with the origin O, where $x$ is appointed vertically along the surface, and $y$ is perpendicular to it. The impact of the heat source/sink and velocity slip is taken into account for the boundary layer and heat transfer analyses. The wall surface is stretched ($b>0$) or shrunk ($b<0$) with surface velocity $u_w\left(x\right)=b{e}^{x/L}$ in the exponential form. $B\left(x\right)$ is a transverse magnetic field that is perpendicular to the stretching/shrinking surface. Here, $T_w\left(x\right)$ is the surface temperature with $T_w\left(x\right)= T_{\infty }+ T_0{e}^{2x/L}$, whereby $T_0<0$ is the constant cooled surface (opposing flow) and $T_0>0$ is the constant heated surface (assisting flow). When the buoyancy forces resulting from the temperature differential between the wall and the free stream increase and considerably alter the flow and thermal fields, the mixed convection in stagnation flows becomes significant, as reported by Ramachandran et al. [24]. This study also makes the crucial assumption that the stagnation point occurs when $T_w\left(x\right)= T_{\infty }$, where the temperature of the plate reaches the same level as the surrounding air. A fluid particle has zero velocity at the stagnation point in the flow field, which is where the fluid region is closest to the stagnation point. The free stream velocity of the flow close to the stagnation point is often given by the equation $U_{\infty }=U_{\infty }\left(x\right)=a{e}^{x/L}$, where $a$ is a constant for $a>0$, $x$ is the distance from stagnation point, and $L$ is denoted as the reference length. The free stream travels in the direction of the plate’s surface before splitting into two to pass it. The velocity at the center of the surface is fully at rest because there is no velocity in either the $x$ or $y$ directions. Here, the center of the surface is referred to as stagnation point.

As for the nanoparticle volume fraction in hybrid ferrofluids (${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4/\mathrm{water}$), ${\mathrm{Fe}}_3{\mathrm{O}}_4$ is the first ferroparticle volume fraction parameter, ${\varphi }_1$, while ${\mathrm{CoFe}}_2{\mathrm{O}}_4$ is the second ferroparticle volume fraction parameter, ${\varphi }_2$. The shape formed for the hybrid ferrofluids is in a stable spherical shape of nanoparticle and uniform size, while the agglomeration effects are disregarded. In this instance, Tiwari and Das’ [25] mathematical nanofluid model is used. It is worth noting that this nanofluid model uses a single-phase approach and assumes that all of the particles are homogeneous in size and shape. Additionally, it ignores interactions between the particles and the fluid they are in contact with [68,69]. Since the base fluid can be thought of as acting as a single fluid when it is readily fluidized, this presumption is realistic and supports the use of the single-phase model in this investigation. The boundary layer is assumed to be laminar, as the fluid flow velocity changes evenly as the fluid flow moves away from the wall surface. Of note, the assumption made on the boundary layer approximation is due to the boundary layer’s edge responding to the external flow in the same way that an object’s actual surface would. Therefore, the boundary layer provides every object with an “effective” shape that is typically slightly distinct from its actual shape. Under the assumption of Boussinesq and boundary layer approximations mentioned above, the governing equations can be written as [22,50].

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,$$

(1)

$$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=U_{\infty }\frac{d{U}_{\infty }}{dx}+\frac{{\mu }_{hnf}}{{\rho }_{hnf}}\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\sigma }_{hnf}{B}^2}{{\rho }_{hnf}}\left(U_{\infty }-u\right)+\frac{{\left(\rho \beta \right)}_{hnf}}{{\rho }_{hnf}}\left(T-T_{\infty }\right)g,$$

(2)

$$u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}={\alpha }_{hnf}\frac{{\partial }^{2}T}{\partial y^2}+\frac{Q_1}{{\left(\rho C_p\right)}_{hnf}}\left(T-T_{\infty }\right),$$

(3)

subject to

$$\begin{array}{c}u=u_w\left(x\right)+{\xi }_1\left(x\right)\frac{\partial u}{\partial y}, v=v_w\left(x\right), T=T_w\left(x\right) \mathrm{at} y=0,\\ u\to U_{\infty }\left(x\right), T\to T_{\infty } \mathrm{as} y\to \infty .\end{array}$$

(4)

At this point, $u$ represents the velocity component along the x-axis, while $v$ represents the velocity component along the y-axis. The acceleration due to gravity is denoted as g, T represents the hybrid ferrofluid temperature, ${\xi }_1\left(x\right)$ refers to the velocity slip factor, and constant heat source/sink is denoted as $Q_1$. Furthermore, ${\left(\rho \beta \right)}_{hnf}$, ${\left(\rho C_p\right)}_{hnf}{\mu }_{hnf}$, ${\rho }_{hnf}$ and ${\sigma }_{hnf}$ portray the thermal expansion coefficient, heat capacitance, dynamic viscosity, density, and electrical conductivity of the hybrid ferrofluids (${\mathrm{Fe}}_3{\mathrm{O}}_4/{\mathrm{CoFe}}_2{\mathrm{O}}_4$), respectively. The thermal diffusivity of the hybrid ferrofluids is described as ${\alpha }_{hnf}=k_{hnf}/{\left(\rho C_p\right)}_{hnf}$, where $k_{hnf}$ denotes the thermal conductivity of the hybrid ferrofluids. Here, the theoretical models for thermophysical properties of ferrofluid and hybrid ferrofluids are presented in

Table 1. Theoretical models for properties of the ferrofluid and hybrid ferrofluids [37,50,70].

Properties	Ferrofluid	Hybrid Ferrofluids
Viscosity	${\mu }_{nf}=\frac{{\mu }_f}{{\left(1-\varphi \right)}^{2.5}}$	${\mu }_{hnf}=\frac{{\mu }_f}{{\left(1-{\varphi }_1\right)}^{2.5}{\left(1-{\varphi }_2\right)}^{2.5}}$
Density	${\rho }_{nf}=\left(1-{\varphi }_1\right){\rho }_f+{\varphi }_1{\rho }_{s1}$	${\rho }_{hnf}=\left(1-{\varphi }_2\right)\left[\left(1-{\varphi }_1\right){\rho }_f+{\varphi }_1{\rho }_{s1}\right]+{\varphi }_2{\rho }_{s2}$
Electrical Conductivity	$\frac{{\sigma }_{nf}}{{\sigma }_f}=1+\frac{3\left({\sigma }_{s1}/{\sigma }_f-1\right){\varphi }_1}{{\sigma }_{s1}/{\sigma }_f+2-\left({\sigma }_{s1}/{\sigma }_f-1\right){\varphi }_1}$	$\frac{{\sigma }_{hnf}}{{\sigma }_{nf}}=\frac{{\sigma }_{s2}+2{\sigma }_{nf}-2{\varphi }_2\left({\sigma }_{nf}-{\sigma }_{s2}\right)}{{\sigma }_{s2}+2{\sigma }_{nf}+{\varphi }_2\left({\sigma }_{nf}-{\sigma }_{s2}\right)}$ where $\frac{{\sigma }_{nf}}{{\sigma }_f}=\frac{{\sigma }_{s1}+2{\sigma }_f-2{\varphi }_1\left({\sigma }_f-{\sigma }_{s1}\right)}{{\sigma }_{s1}+2{\sigma }_f+{\varphi }_1\left({\sigma }_f-{\sigma }_{s1}\right)}$
Thermal Diffusivity	${\alpha }_{nf}=\frac{k_{nf}}{{\left(\rho C_p\right)}_{nf}}$	${\alpha }_{hnf}=\frac{k_{hnf}}{{\left(\rho C_p\right)}_{hnf}}$
Thermal Expansion Coefficient	${\left(\rho \beta \right)}_{nf}=\left(1-{\varphi }_1\right){\left(\rho \beta \right)}_f+{\varphi }_1{\left(\rho \beta \right)}_{s1}$	${\left(\rho \beta \right)}_{hnf}=\left(1-{\varphi }_2\right)\left[\left(1-{\varphi }_1\right){\left(\rho \beta \right)}_f+{\varphi }_1{\left(\rho \beta \right)}_{s1}\right]+{\varphi }_2{\left(\rho \beta \right)}_{s2}$
Heat Capacity	${\left(\rho C_p\right)}_{nf}=\left(1-{\varphi }_1\right){\left(\rho C_p\right)}_f+{\varphi }_1{\left(\rho C_p\right)}_{s1}$	${\left(\rho C_p\right)}_{hnf}=\left(1-{\varphi }_2\right)\left[\left(1-{\varphi }_1\right){\left(\rho C_p\right)}_f+{\varphi }_1{\left(\rho C_p\right)}_{s1}\right]+{\varphi }_2{\left(\rho C_p\right)}_{s2}$

, and

Table 2. Thermophysical properties of ferroparticles and water (base fluid) [13,18].

Thermophysical Properties	$F{e}_3{O}_4$	$CoF{e}_2{O}_4$	$H_{2}O$
$C_p\left(Jk{g}^{-1}{K}^{-1}\right)$	670	700	4179
$\rho \left(kg{m}^{-3}\right)$	5180	4907	997.1
$k\left(W{m}^{-1}{K}^{-1}\right)$	9.7	3.7	0.613
$\beta \times {10}^{-5}\left(K^{-1}\right)$	1.3	1.3	21
$\sigma \left(S{m}^{-1}\right)$	$0.74\times {10}^6$	$1.1\times {10}^7$	$5.5\times {10}^{-6}$

shows the physical properties of nanoparticles ${\mathrm{Fe}}_3{\mathrm{O}}_4$ and ${\mathrm{CoFe}}_2{\mathrm{O}}_4$ together with ${\mathrm{H}}_2\mathrm{O}$ particles.

For obtaining the similarity solutions of the governing Equations (1)–(3) subject to Equation (4), the terms $v_w\left(x\right)$, $B\left(x\right)$, $Q_1\left(x\right)$ and ${\xi }_1\left(x\right)$ are given in the following form [26,50,52]:

$$v_w\left(x\right)= -\sqrt{\frac{a{v}_f}{2L}}{e}^{x/2L}s, B\left(x\right)=B_0{e}^{x/2L}, Q_1\left(x\right)=Q_0{e}^{x/L}, {\xi }_1\left(x\right)={\xi }_0{e}^{-x/2L},$$

(5)

where continuous mass flux velocity is denoted as $s$ with suction condition $(s>0)$, injection condition $(s<0)$, and impenetrable surface $\left(s=0\right),$ while $B_0$, $Q_0$ and ${\xi }_0$ are constants. Next, an appropriate similarity transformation is considered as follows by pursuing [50].

$$\begin{array}{c}u=a{e}^{x/L}{f}^{\prime }\left(\eta \right), v=-\sqrt{\frac{a{v}_f}{2L}}{e}^{x/2L} \left(f\left(\eta \right)+\eta f^{\prime }\left(\eta \right)\right),\\ \theta \left(\eta \right)=\frac{T-T_{\infty }}{ T_w-T_{\infty }}, \eta =\sqrt{\frac{a}{2v_{f}L}}{e}^{x/2L}y,\end{array}$$

(6)

where $f^{\prime }$ indicates differentiation concerning $\eta$. Here, the reduction of the differential equation will be made by implying Equations (5) and (6) into Equations (2) and (3) along with Equation (4). A reduced form of governing boundary layer system is obtained as follows:

$$\frac{{\mu }_{hnf}}{{\mu }_f}{f}^{\prime \prime \prime }+\frac{{\rho }_{hnf}}{{\rho }_f}\left(f{f}^{″}-2{\left(f^{\prime }\right)}^2+2\right)+\frac{{\sigma }_{hnf}}{{\sigma }_f}M\left(1-f^{\prime }\right)+2\frac{{\left(\rho \beta \right)}_{hnf}}{{\left(\rho \beta \right)}_f}\lambda \theta =0,$$

(7)

$$\frac{1}{\mathrm{Pr}}\frac{k_{hnf}}{k_f}{\theta }^{″}+\frac{{\left(\rho C_p\right)}_{hnf}}{{\left(\rho C_p\right)}_f}\left(f{\theta }^{\prime }-4f^{\prime }\theta \right)+Q\theta =0,$$

(8)

subject to

$$\begin{array}{c}f\left(\eta \right)=s, f^{\prime }\left(\eta \right)=c+\xi f^{″}\left(\eta \right), \theta \left(\eta \right)=1 \mathrm{at} \eta =0\\ f^{\prime }\left(\eta \right)=1, \theta \left(\eta \right)=0 \mathrm{as} \eta \to \infty .\end{array}$$

(9)

Here, $\lambda$ represents the mixed convection parameter with the case opposing flow $(\lambda <0)$ and assisting flow $(\lambda >0)$, $\mathrm{Pr}$ indicates the Prandtl number, $M$ denotes the magnetic parameter, $Q$ represents the presence of heat source $(Q>0)$ and heat sink $(Q<0)$, $c$ denotes stretching $(c>0)$ and shrinking $(c<0)$ parameter on the surface, and $\xi$ is the velocity slip parameter, which is defined as follows:

$$M=\frac{2L{B}_0^2{\sigma }_f}{a{\rho }_f}, \lambda =\frac{G{r}_x}{{\mathrm{Re}}_x^2}, \mathrm{Pr}=\frac{{\left(\mu C_p\right)}_f}{k_f}, Q=\frac{2L{Q}_0}{a{\left(\rho C_p\right)}_f}, c=\frac{b}{a}, \xi ={\xi }_0\sqrt{\frac{a}{2v_{f}L}},$$

(10)

where the local Reynolds number ${\mathrm{Re}}_x$ and the local Grashof number $G{r}_x$ are given by

$$G{r}_x=\frac{g{\beta }_f\left(T_w-T_{\infty }\right)L^3}{v_f}, {\mathrm{Re}}_x=\frac{L{U}_{\infty }\left(x\right)}{{\nu }_f}.$$

(11)

Next, the skin friction coefficient $C_f$ and the local Nusselt number $N{u}_x$, as the physical quantities of practical interest in this problem, are defined as follows [50]:

$$C_f=\frac{{\tau }_w}{{\rho }_f{U}_{\infty }^2}, N{u}_x=\frac{L{q}_w}{k_f\left(T_w-T_{\infty }\right)},$$

(12)

where $q_w$ is the surface heat flux and ${\tau }_w$ is the wall shear stress, as defined below:

$${\tau }_w={\mu }_{hnf}{\left(\frac{\partial u}{\partial y}\right)}_{y=0}, q_w=-k_{hnf}{\left(\frac{\partial T}{\partial y}\right)}_{y=0}.$$

(13)

Using Equation (6) in Equations (12) and (13), the following is obtained:

$$\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f=\frac{{\mu }_{hnf}}{{\mu }_f}{f}^{″}\left(0\right), \sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x=-\frac{k_{hnf}}{k_f}{\theta }^{\prime }\left(0\right).$$

(14)

3. Stability Analysis

There is an appearance of dual solutions from the numerical outputs of Equations (7) and (8) with respective governing parameters. Hence, a stability analysis was conducted to identify a stable solution (first solution) and an unstable solution (second solution). The unsteady forms of Equations (1)–(4) are assumed by introducing a new similarity transformation cooperatively with the similarity solutions (Equation (5)) for the unsteady problem:

$$\begin{array}{c}u=a{e}^{x/L}\frac{\partial f}{\partial \eta }\left(\eta ,\tau \right), v=-\sqrt{\frac{a{\nu }_f}{2L}}{e}^{x/2L}\left[f\left(\eta ,\tau \right)+\eta \frac{\partial f}{\partial \eta }+2\tau \frac{\partial f}{\partial \tau }\right],\\ \theta \left(\eta ,\tau \right)= \frac{T-T_{\infty }}{ T_w-T_{\infty }}, \eta =\sqrt{\frac{a}{2v_{f}L}}{e}^{x/2L}y, \tau =\frac{a}{2L}t{e}^{x/L}\end{array}$$

(15)

Here, Equation (1) was held, and Equations (2) and (3) were formed as:

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=U_{\infty }\frac{d{U}_{\infty }}{dx}+ \frac{{\mu }_{hnf}}{{\rho }_{hnf}}\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\sigma }_{hnf}{B}^2}{{\rho }_{hnf}}\left(U_{\infty }-u\right)+\frac{{\left(\rho \beta \right)}_{hnf}}{{\rho }_{hnf}}\left(T-T_{\infty }\right)g,$$

(16)

$$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}={\alpha }_{hnf}\frac{{\partial }^{2}T}{\partial y^2}+\frac{Q_1}{{\left(\rho C_p\right)}_{hnf}}\left(T-T_{\infty }\right),$$

(17)

where t denotes the time. Thus, the transformed Equations (16) and (17) may be written as follows:

$$\begin{array}{c}\frac{{\mu }_{hnf}}{{\mu }_f}\frac{{\partial }^{3}f}{\partial {\eta }^3}+\frac{{\rho }_{hnf}}{{\rho }_f}\left[f\frac{{\partial }^{2}f}{\partial {\eta }^2}-2{\left(\frac{\partial f}{\partial \eta }\right)}^2+2-\frac{{\partial }^{2}f}{\partial \tau \partial \eta }-2\tau \left[\frac{{\partial }^{2}f}{\partial \tau \partial \eta }\frac{\partial f}{\partial \eta }-\frac{{\partial }^{2}f}{\partial {\eta }^2}\frac{\partial f}{\partial \tau }\right]\right]+\frac{{\sigma }_{hnf}}{{\sigma }_f}M\left[1-\frac{\partial f}{\partial \eta }\right]\\ +2\lambda \frac{{\left(\rho \beta \right)}_{hnf}}{{\left(\rho \beta \right)}_f}\theta =0,\end{array}$$

(18)

$$\frac{1}{Pr}\frac{k_{hnf}}{k_f}\frac{{\partial }^2\theta }{\partial {\eta }^2}+2Q\theta +\frac{{\left(\rho C_p\right)}_{hnf}}{{\left(\rho C_p\right)}_f}\left[f\frac{\partial \theta }{\partial \eta }-4\theta \frac{\partial f}{\partial \eta }-\frac{\partial \theta }{\partial \tau }-2\tau \left[\frac{\partial f}{\partial \eta }\frac{\partial \theta }{\partial \tau }-\frac{\partial f}{\partial \tau }\frac{\partial \theta }{\partial \eta }\right]\right]=0,$$

(19)

subject to

$$\begin{array}{c}f\left(\eta ,\tau \right)+2\tau \frac{\partial f}{\partial \tau }=s, \frac{\partial f}{\partial \eta }\left(\eta ,\tau \right)=c+\xi \frac{{\partial }^{2}f}{\partial {\eta }^2}\left(\eta ,\tau \right), \theta \left(\eta ,\tau \right)=1 \mathrm{at} \eta =0,\\ \frac{\partial f}{\partial \eta }\left(\eta ,\tau \right)=1, \theta \left(\eta ,\tau \right)=0 \mathrm{as} \eta \to \infty .\end{array}$$

(20)

Following Weidman and Turner for analyzing the stability of dual solutions [67], we write

$$f\left(\eta ,\tau \right)=f_0\left(\eta \right)+e^{-\gamma \tau }F\left(\eta ,\tau \right), \theta \left(\eta ,\tau \right)={\theta }_0\left(\eta \right)+e^{-\gamma \tau }G\left(\eta ,\tau \right),$$

(21)

where $\gamma$ denotes the unknown eigenvalue that decides the stability of the solutions, whereas $F\left(\eta \right)$ and $G\left(\eta \right)$ are comparatively small compared to $f_0\left(\eta \right)$ and ${\theta }_0\left(\eta \right)$. The exponential disturbance demonstrates the rapid decline or the development of the disturbance. Hence, the following expressions can be obtained by implementing time-dependent solutions of Equation (21):

$$\begin{array}{c}\frac{{\mu }_{hnf}}{{\mu }_f}\frac{{\partial }^{3}F}{\partial {\eta }^3}+\frac{{\rho }_{hnf}}{{\rho }_f}[f_0\frac{{\partial }^{2}F}{\partial {\eta }^2}+f_0^{\prime \prime }F-4f_0^{\prime }\frac{\partial F}{\partial \eta }+\gamma \frac{\partial F}{\partial \eta }-\frac{{\partial }^{2}F}{\partial \eta \partial \tau }\\ -2\tau \left[-\gamma f_0^{\prime }\frac{\partial F}{\partial \eta }+f_0^{\prime }\frac{{\partial }^{2}F}{\partial \eta \partial \tau }+\gamma f_0^{\prime \prime }F-f_0^{\prime \prime }\frac{\partial F}{\partial \tau }\right]]-\frac{{\sigma }_{hnf}}{{\sigma }_f}M\frac{\partial F}{\partial \eta }+2\lambda \frac{{\left(\rho \beta \right)}_{hnf}}{{\left(\rho \beta \right)}_f}G=0,\end{array}$$

(22)

$$\begin{array}{ll}\frac{1}{Pr}\frac{k_{hnf}}{k_f}\frac{{\partial }^{2}G}{\partial {\eta }^2}& +2QG\\ & +\frac{{\left(\rho C_p\right)}_{hnf}}{{\left(\rho C_p\right)}_f}[\left(\frac{\partial G}{\partial \eta }{f}_0+F{\theta }_0^{\prime }-4\frac{\partial F}{\partial \eta }{\theta }_0-4G{f}_0^{\prime }+\gamma G-\frac{\partial G}{\partial \tau }\right)\\ & -2\tau \left[\left(-f_0^{\prime }\gamma G+\frac{\partial G}{\partial \tau }{f}_0^{\prime }+\gamma F{\theta }_0^{\prime }-\frac{\partial F}{\partial \tau }{\theta }_0^{\prime }\right)\right]]=0,\end{array}$$

(23)

subject to

$$\begin{array}{c}F\left(\eta ,\tau \right)+\frac{2\tau }{1-2\tau \gamma }\frac{\partial F}{\partial \tau }\left(\eta ,\tau \right)=0, \frac{\partial F}{\partial \eta }\left(\eta ,\tau \right)=\xi \frac{{\partial }^{2}F}{\partial {\eta }^2}\left(\eta ,\tau \right), G\left(\eta ,\tau \right)=0 \mathrm{at} \eta =0,\\ \frac{\partial F}{\partial \eta }=0, G\left(\eta ,\tau \right)=0 \mathrm{as} \eta \to \infty .\end{array}$$

(24)

To point out the early growth/decay of the solution, in Equation (21), $\tau$ is set at 0, which yields $F\left(\eta ,\tau \right)=F_0\left(\eta \right)$ and $G\left(\eta ,\tau \right)=G_0\left(\eta \right)$. Hence, via the collocation method, the following linearized eigenvalue problem was resolved:

$$\frac{{\mu }_{hnf}}{{\mu }_f}{F}_0^{\prime \prime \prime }+\frac{{\rho }_{hnf}}{{\rho }_f}\left[f_0{F}_0^{\prime \prime }+f_0^{\prime \prime }{F}_0-4f_0^{\prime }{F}_0^{\prime }+\gamma F_0^{\prime }\right]-\frac{{\sigma }_{hnf}}{{\sigma }_f}M{F}_0^{\prime }+2\lambda \frac{{\left(\rho \beta \right)}_{hnf}}{{\left(\rho \beta \right)}_f}{G}_0=0,$$

(25)

$$\frac{1}{Pr}\frac{k_{hnf}}{k_f}{G}_0^{\prime \prime }+\frac{{\left(\rho C_p\right)}_{hnf}}{{\left(\rho C_p\right)}_f}\left[G_0^{\prime }{f}_0+F_0{\theta }_0^{\prime }-4G_0{f}_0^{\prime }-4F_0^{\prime }{\theta }_0+\gamma G_0\right]+2Q{G}_0=0,$$

(26)

accompanied by

$$\begin{array}{c}{F}_0\left(\eta \right)=0, F{\prime }_0\left(\eta \right)=\xi F_0^{\prime \prime }\left(\eta \right), G_0\left(\eta \right)=0 \mathrm{at} \eta =0,\\ F_0^{\prime }\left(\eta \right)=0, G_0\left(\eta \right)=0 \mathrm{as} \eta \to \infty .\end{array}$$

(27)

The infinite set of eigenvalues ${\gamma }_1<{\gamma }_2<{\gamma }_3<$… is provided by Equations (25)–(27). In determining the stability of the dual solution, it is important to obtain the smallest eigenvalue, ${\gamma }_1$. As reported by Weidman and Turner [67], one of the boundary conditions $F_0^{\prime }(\infty )=0$ or $G_0(\infty )=0$ will be relaxed, displaying the infinite set of eigenvalues ${\gamma }_1<{\gamma }_2<{\gamma }_3<$…. Here, the boundary condition $F_0^{\prime }(\infty )=0$ was selected to be relaxed and substituted with a new boundary condition $F_0^{\prime \prime }\left(0\right)=1$ to solve the linearized equation mentioned above.

4. Results and Discussion

In evaluating Equations (7) and (8), the numerical computation employed the bvp4c solver in MATLAB software to obtain the numerical solution for the formulated boundary value problem along with the boundary conditions in Equation (9). In this present work, the influence of various governing parameters, such as nanoparticle volume fraction ${\varphi }_1$, ${\varphi }_2$, magnetic parameter M, heat source/sink parameter Q, mixed convection parameter $\lambda$, velocity slip $\xi$ and stretching/shrinking surface parameter c, are considered in the problem. During numerical computation in MATLAB, values ${\varphi }_1=0.01$, ${\varphi }_2=0.01$, M = 0.1, Q = 0.1 $\xi =1$ and c = −0.5 are fixed and can be changed based on various governing parameters selected. For validating the present model, the results were compared with Waini et al. [50], Bachok et al. [71], and Ur Rehman et al. [72] for nanofluids, ${\mathrm{Al}}_2{\mathrm{O}}_3$/water between skin friction coefficient, $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ and Nusselt number, $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x.$ These were computed as in

Table 3. Comparison values of $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ and $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ for ${\mathrm{Al}}_2{\mathrm{O}}_3$/water (${\varphi }_2=0$) with various values of c and ${\varphi }_1$ when $\lambda =0$, Pr = 6.2.

		$$\sqrt{2} R{e}_x^{1/2}{C}_f$$				$$\sqrt{2} R{e}_x^{-1/2}N{u}_x$$
c	$${\varphi }_1$$		Rehman et al. [72]	Waini et al. [50]	Present Results	Waini et al. [50]	Present Results
−0.5	0	2.1182		2.1182	2.118169	0.0588	0.058787
	0.1	2.7531		2.7531	2.753091	0.4439	0.443904
	0.2	3.5372		3.5372	3.537175	0.7636	0.763593
0	0	1.6872	1.68720	1.6872	1.687218	2.5066	2.506625
	0.1	2.1929	2.19293	2.1930	2.192963	2.9655	2.965516
	0.2	2.8174	2.81750	2.8175	2.817522	3.4292	3.429219
0.5	0	0.9604	0.96040	0.9604	0.960416	4.0816	4.081573
	0.1	1.2483	1.24829	1.2483	1.248302	4.6637	4.663651
	0.2	1.6039	1.60399	1.6038	1.603819	5.2726	5.272615

, which shows positive comparative results.

Figure 2 and Figure 3 display numerical results for the effect of several values of nanoparticle volume fraction ${\varphi }_2$ on variations of $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ and $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ against $\lambda$. Both figures illustrate that the values $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ and $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ for hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0.02$) are the highest compared to hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0.01$) and ferrofluid, ${\mathrm{Fe}}_3{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0$). Physically, the presence of ${\varphi }_2$ led to a collision between two solid ferroparticles scattered in the base fluid, which led to a greater value of $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ than for one solid ferroparticle. The presence of ${\mathrm{CoFe}}_2{\mathrm{O}}_4$ as an additional ferroparticle does influence the performance of the heat transfer rate between hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water and ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$/water. In terms of heat transfer, as in Figure 3, hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$- ${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water is better than ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$/water. It is worth noticing that the increment in ${\varphi }_2$ values also contributes to the enhancement of $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ and $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$, as portrayed in Figure 2 and Figure 3, respectively. Increasing the nanoparticle volume fraction will enhance its synergistic effect, thereby enhancing the heat transfer rate. From the observation, for assisting flow region ($\lambda >0$), the values of $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ and $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ are greater than in the opposing flow region ($\lambda <0$). Interestingly, the existence of a dual solution can be observed in the range $\lambda >{\lambda }_c$, where ${\lambda }_c$ denotes the critical value of $\lambda$. Note that the boundary layer separation occurs between the first and second solutions at $\lambda ={\lambda }_c$ (the point of bifurcation). It can be observed in Figure 2 and Figure 3 that the ${\lambda }_c$ for hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water is lower than that for ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$/water. This implies that hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water can delay the boundary layer separation and minimize the range of solutions compared with ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$/water.

Figure 4 delineates the magnetic parameter $M$ impact on $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ with $\lambda$. The increased value of $M$ demonstrates the increased value for $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$. This phenomenon occurs because the presence and upsurge value of $M$ increases the Lorentz power, otherwise called resistive sort power. This creates significantly more resistance from the flow and increases the shear stress at the wall. Other observations show the presence of $M$ decreased the val of ${\lambda }_c,$. which likely slows down the separation of the boundary layer. The variation of $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ against $\lambda$ with selected values of heat source Q (> 0) and heat sink Q (< 0) are portrayed in Figure 5. For the increment strength of Q (> 0), $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ showed a decrease in value, while the increment strength of Q (< 0) caused the values of $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ to increase. The observation on Q (> 0) indicates a reduction in heat flux because there is an increased thickness of the thermal boundary layer when there is an additional heat source in the boundary layer. Meanwhile, the heat flux was enhanced with the increment of Q (< 0) because of the reduced thickness of the thermal boundary layer. It was believed that as the strength of Q (> 0) increased, the value of ${\lambda }_c$ reduced, while as the strength of Q (< 0) increased, the value of ${\lambda }_c$ rose, as seen in Figure 5. This observation suggests that as the strength of Q (> 0) grows, it may slow down the boundary layer separation, while as the strength of Q (< 0) grows, it may accelerate the boundary layer separation.

The plots of $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ against the mixed convection parameter $\lambda$ with numerous values of velocity slip parameter $\xi$ are visualized in Figure 6. The trend of $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ decreases with the escalation of $\xi$ value for the first solution in both assisting ($\lambda >0$) and opposing flow ($\lambda <0$). In the physical aspect, the velocity slip occurs due to the finite number of molecular interactions between particles from the fluid and the wall. Generally, the fluid flow velocity is influenced by the presence of velocity slip. Here, the presence of the velocity slip parameter can be used as a controlling factor to obtain the outcome needed, whether to have a better or slower fluid flow velocity. For observation of boundary layer separation, the ${\lambda }_c$ decreased in value with the increase in $\xi$, allowing the conclusion that an increment of $\xi$ may delay the boundary layer separation.

Table 4. $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ values for various values of c.

$$\sqrt{2} R{e}_x^{1/2}{C}_f$$
λ	c
	−1.5	−1.0	−0.5	0.5	1.0	1.5
5	2.168550	1.784748	1.397650	0.614603	0.219096	−0.178813
4	2.115534	1.734606	1.350040	0.571247	0.177546	−0.218730
3	2.060326	1.682616	1.300853	0.526709	0.134958	−0.259564
2	2.002605	1.628538	1.249907	0.480876	0.091243	−0.301388
1	1.941959	1.572076	1.196980	0.433617	0.046297	−0.344285
0	1.877853	1.512850	1.141798	0.384778	0.000000	−0.388351
−1	1.809554	1.450369	1.084017	0.334170	−0.047792	−0.433696
−2	1.736025	1.383972	1.023190	0.281566	−0.097249	−0.480452
−3	1.655703	1.312731	0.958725	0.226683	−0.148578	−0.528774
−4	1.566017	1.235271	0.889801	0.169164	−0.202032	−0.578847
−5	1.462139	1.149388	0.815224	0.108540	−0.257932	−0.630900

and

Table 5. $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ values for various values of c.

$$\sqrt{2} R{e}_x^{-1/2}N{u}_x$$
λ	c
	−1.5	−1.0	−0.5	0.5	1.0	1.5
5	8.570820	8.855438	9.123041	9.615376	9.843333	10.060869
4	8.441423	8.738077	9.015792	9.524113	9.758478	9.981619
3	8.303596	8.613944	8.903014	9.429008	9.670365	9.899579
2	8.155904	8.482027	8.783971	9.329647	9.578673	9.814497
1	7.996487	8.341045	8.657756	9.225537	9.483029	9.726082
0	7.822852	8.189350	8.523232	9.116087	9.382989	9.633998
−1	7.631547	8.024754	8.378948	9.000581	9.278027	9.537849
−2	7.417554	7.844264	8.222997	8.878132	9.167509	9.437166
−3	7.173119	7.643590	..052800	8.747629	9.050654	9.331386
−4	6.885144	7.416234	7.864719	8.607640	8.926494	9.219827
−5	6.528230	7.151525	7.653359	8.456280	8.793795	9.101644

present the $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ and $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ against $\lambda$ for several values of shrinking parameter ($c<0$) and stretching parameter ($c>0$), respectively. The values of $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ and $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ in assisting flow region ($\lambda >0$) are higher than in opposing flow region ($\lambda <0$) for both stretching/shrinking parameters. Additionally, the decreased value of the shrinking parameter ($c<0$) gives a better $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ but reduces $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$. An observation suggests that the velocity flow will increase with the shear stress at the wall, while a reduction in heat transfer tends to occur. Meanwhile, the increased value of the stretching parameter ($c>0$) shows a decreasing trend value for $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ but an increase in $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$. Here, the velocity flow tends to slow down and reduce the shear stress, while the increment of heat transfer tends to occur.

The profiles of velocity $f\prime \left(\eta \right)$ and temperature $\theta \left(\eta \right)$ in Figure 7 and Figure 8, respectively, highlight the influence of nanoparticle volume fraction ${\varphi }_2$ (for ${\mathrm{CoFe}}_2{\mathrm{O}}_4$) in opposing flow region $\left(\lambda =-2\right)$ with ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0$), hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0.01$) and hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0.02$). Interestingly, it is depicted in Figure 7 that hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0.02$) has a slightly greater $f^{\prime }\left(\eta \right)$ than ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0$) and hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0.01$). It is clear that the presence of ${\varphi }_2$ results in an increase in the $f\prime \left(\eta \right)$ but a reduction in momentum boundary layer thickness. Moreover, the increment in values of ${\varphi }_2$ yields an enhancement in $f\prime \left(\eta \right)$. Physically, this will increase skin friction since the shear stress is enhanced because of the reduction in the momentum of the boundary layer thickness. This outcome is consistent with Figure 2, which shows that hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0.02$) has a higher $\sqrt{2}{ \mathrm{Re}}_x^{1/2}{C}_f$ compared to others in opposing flow regions $\left(\lambda =-2\right)$. Regarding the behavior of $\theta \left(\eta \right)$ as shown in Figure 8, ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0$) can be seen to have a slightly lower $\theta \left(\eta \right)$ than hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0.01$) and hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0.02)$. This discovery leads to the conclusion that the thermal boundary layer thickness concerning ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0$) is marginally less compared to that of the hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water (${\varphi }_1=0.01$, ${\varphi }_2=0.01$) and hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4/\mathrm{water}$ (${\varphi }_1=0.01$, ${\varphi }_2=0.02$. Physically, thinner thermal boundary layer thickness led to a larger heat flux.

The behaviors of $f^{\prime }\left(\eta \right)$ and $\theta \left(\eta \right)$ with the adjustment of the heat source Q (> 0) and heat sink Q (< 0) in opposing flow region $\left(\lambda =-2\right)$ are shown in Figure 9 and Figure 10, respectively. As portrayed in Figure 9, the increment in strength of heat source Q (> 0) slightly decreases the velocity profile, which in turn increases the boundary layer thickness. It slows down the fluid flow velocity, which in turn reduces the surface shear stress. As for the increment in strength of heat sink Q (< 0), it gives the converse result, which increased the velocity and shear stress while reducing boundary layer thickness. Meanwhile, the value of $\theta \left(\eta \right)$ and thermal boundary layer thickness as displayed in Figure 10 increase with the increment in strength of heat source Q (> 0), while the increased strength of heat sink Q (< 0) shows the reverse result. Here, the increment in strength of heat source Q (> 0) physically reduces the heat flux on the surface with the built-up thermal boundary layer thickness. Meanwhile, there exists a contrary statement for the increased strength of heat sink Q (< 0), where the heat flux on the surface is physically increased due to the reduction in the thermal boundary layer thickness. This result aligns with Figure 5, where the value of $\sqrt{2}{ \mathrm{Re}}_x^{-1/2}N{u}_x$ decreases with the increment in strength of heat source Q (> 0) while it increases with the increment in strength of heat sink Q (< 0).

The existence of a non-unique solution (dual solution) in the assisting ($\lambda >0$) and opposing flow ($\lambda <0$) region from Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 can be observed. Generally, in the range of $\lambda >{\lambda }_c$, the dual solution will appear. A unique solution where the first and second solution meet is known as the critical value ${\lambda }_c$. Point $\lambda ={\lambda }_c$ is known as a point of bifurcation of first and second solutions. Separation for the boundary layer occurs when $\lambda ={\lambda }_c$, and using boundary layer approximation, no solution is obtained if the range is between $\lambda <{\lambda }_c$. Therefore, it is necessary to perform the flow’s stability analysis. In order to execute the stability analysis, the same numerical approach was used for solving the linearized equations (25)–(27) numerically, as discussed in the previous section. It is essential to obtain small eigenvalues $\gamma$ for verification of stability analysis. Here,

Table 6. The smallest eigenvalue $\gamma$ with the variation of $\lambda$ when Pr = 6.2, ${\varphi }_1=0.01$, ${\varphi }_2=0.01$, $c=-0.5$, $\xi =1$ and $s=1$ for ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water.

$$Q$$	$$M$$	$$\lambda$$	$${\gamma }_1$$	$${\gamma }_2$$
0.1	0	−6 −9.1 −9.43	3.6558 1.2011 0.4997	−3.3488 −1.2430 −0.5765
0.1	0.5	−9.4 −10.4 −10.52	2.0449 0.6352 0.0656	−2.0046 −0.7137 −0.1580
0.1	1.0	−9 −10.4 −11.44	3.1228 2.0376 0.4799	−2.8976 −1.9977 −0.5702

was constructed to demonstrate the behavior of small eigenvalues $\gamma$ when M = 0.1, Q = 1, c = −0.5, $\xi$ = 0.1, and of hybrid ferrofluids (${\varphi }_1$ = 0.01, ${\varphi }_2$ = 0.01) concerning $\lambda$. According to [66] and [67], the first solution gives a positive value for eigenvalue ${\gamma }_1,$ which implies the initial decay of disturbance for real and stable flow as time passes. Meanwhile, the negative value of the eigenvalue ${\gamma }_2$ from the second solution implies initial growth of disturbance for unreal and unstable flow as time passes. It can be seen from Table 6 that as mixed convection parameter $\lambda$ approaches the critical values ${\lambda }_c$, the values of small eigenvalues ${\gamma }_1$ and ${\gamma }_2$ approach zero. This allows the conclusion that the first solution is in a stable state, while the second solution is in an unstable state.

5. Conclusion

This study investigated the steady mixed convection in hybrid ferrofluid flow past an exponentially shrinking/stretching surface in a stagnation-point area. The governing equations were solved numerically by applying MATLAB’s built-in solver (bvp4c). The influence of variation of governing parameters such as nanoparticle volume fraction, velocity slip parameter, magnetic parameter, heat source/sink parameter, and stretching/shrinking parameter were discussed in detail. Here, the significant outcomes can be summarized as follows:

• The appearance of a non-unique solution or dual solution for a certain governing parameter in both assisting flow and opposing flow. The bifurcation point occurs at the opposing flow region.
• The stability analysis results show a stable first solution and an unstable second solution.
• The increment of nanoparticle volume fraction ${\varphi }_2$ (${\mathrm{CoFe}}_2{\mathrm{O}}_4$) increases heat transfer and skin friction.
• The addition of ${\mathrm{CoFe}}_2{\mathrm{O}}_4$ nanoparticle to hybrid ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$-${\mathrm{CoFe}}_2{\mathrm{O}}_4$/water enhances the heat transfer rate better than ferrofluid ${\mathrm{Fe}}_3{\mathrm{O}}_4$/water due to the collision between two nanoparticles enhancing the synergistic effects.
• The influence of the magnetic field $M$ enhanced the fluid flow velocity due to the presence of Lorentz forces in flow and slowing down the boundary layer separation.
• The increased strength of the heat source $\left(Q>0\right)$ reduced the heat transfer efficiency because of the larger thermal boundary layer thickness due to the additional external heat sources. Meanwhile, the increased strength of the heat sink $\left(Q<0\right)$ yields the reverse result, which in turn reduced thermal boundary layer thickness due to the heat loss for the heat sink.
• The presence of velocity slip $\xi$ reduced the flow velocity and delayed the boundary layer separation.
• The increment of stretching parameter $\left(c>0\right)$ reduced the skin friction value, which slows down the fluid flow velocity. Meanwhile, the increment of shrinking parameter $\left(c<0\right)$ increases the skin friction value, which accelerates the fluid flow velocity.
• In terms of heat transfer, the increment of stretching parameter $\left(c>0\right)$ gives a better value of the Nusselt number, while the increment of shrinking parameter $\left(c<0\right)$ gives a reverse result.