Rotating 3D Flow of Hybrid Nanofluid on Exponentially Shrinking Sheet: Symmetrical Solution and Duality

Abstract

This article aims to study numerically the rotating, steady, and three-dimensional (3D) flow of a hybrid nanofluid over an exponentially shrinking sheet with the suction effect. We considered water as base fluid and alumina ($A{l}_2{O}_3$), and copper ($Cu$) as solid nanoparticles. The system of governing partial differential equations (PDEs) was transformed by an exponential similarity variable into the equivalent system of ordinary differential equations (ODEs). By applying a three-stage Labatto III-A method that is available in bvp4c solver in the Matlab software, the resultant system of ODEs was solved numerically. In the case of the hybrid nanofluid, the heat transfer rate improves relative to the viscous fluid and regular nanofluid. Two branches were obtained in certain ranges of the involved parameters. The results of the stability analysis revealed that the upper branch is stable. Moreover, the results also indicated that the equations of the hybrid nanofluid have a symmetrical solution for different values of the rotation parameter ($\mathsf{\Omega }$).

1. Introduction

For a number of industrial uses, such as the manufacture of rubber pads, the flow of fluid on the shrinking sheet must be considered. In the production phase, the moving surface is supposed to be compressed to its plane and the shrinking sheet connects with the surrounding fluid both thermally and mechanically. The behavior of a shrinking surface may occur in several materials with specific strengths. Initially, Sakiadis [1] proposed the idea of a flow of the boundary layer on a stretching sheet. Later, Crane [2] modified the concept of Sakiadis and applied it to both exponential and linear stretching surfaces. Recently, the flow over a stretching sheet has received a lot of consideration. Several recent studies [3,4,5,6,7,8] have been conducted in this respect, in which numerous impacts have been examined. Due to the high demand and applications of the shrinking surface, we have considered a 3D flow on the shrinking sheet in this research.

Several researchers have been focusing on nanofluid analysis due to the problem of improving the rate of heat transfer. This fluid can be characterized as a homogenous mixture of nanoparticles and conventional fluids. This concept of a mixture of solid particles and fluid was introduced by Maxwell in the 19th century as an attempt to improve the thermal conductivity of fluids. Compared to the micro-particles, nanoparticles remain dispersed longer and continue in dispersion almost indefinitely if they are below the threshold level or enhanced with the surface. It is observed that various approaches to nanofluid analysis have been adapted, but the most efficient approaches are computational and experimental approaches.

However, the experimental approach is very costly, and therefore the computational approach is preferred to investigate nanofluid by using various models. Initially, the results of these models were compared to the experimental results and found to be in excellent agreement. Hassan et al. [9] used copper oxide particles and concluded that “when particles are added in fluid, convection heat transfer rate is improved but flow velocity is declined”. Naramgari and Sulochana [10] considered nanofluid on the exponential sheet and obtained two solutions by using a Buongiorno model of the nanofluid. Tiwari and Das’ model of a single-phase nanofluid was considered by Dero et al. [11] during the examination of the Casson based nanofluid. They concluded that "velocity profiles and corresponding boundary layer thicknesses decrease by a suspension of nanoparticles of silver and copper, whereas the silver nanoparticles show the greater rate of heat transfer enhancement as compared to copper nanoparticles when suspended in Casson fluid”. Biswakarma et al. [12] examined the aluminum oxide water-based nanofluid and concluded that the heat transfer coefficient is enhanced by as much as 13.8% with the nanofluid. Further, Giri et al. [13] considered the fluid flow in the vertical channel in which they found that both the Nusselt number for local sensitive heat and the Nusselt number for local condensation decline monotonically along the axial direction. Later, Giri et al. [14] used the simpler algorithm to investigate the governing equations of the model. Some interesting outcomes of the nanofluid for various physical conditions and effects can be found in these papers [15,16,17].

As mentioned earlier, due to the increase in the heat transfer rate, many academics and scholars are interested in studying nanofluid. Nanofluid has many applications in new technology eras such as bio-labeling, biocatalysts, biosensors, transportation, biomolecules separation and purification, engine cooling, vehicle thermal management, thermal storage, cooling in nuclear systems, solar water heating, production of glass fiber, defense, and drug delivery. Due to the rising demand for the heat transfer rate from various sectors of the industries, researchers have been attempting to mix many solid nanoparticles with various kinds of base fluids which leads to the discovery of a “hybrid nanofluid” as the new kind of nanofluid. Waini et al. [18] examined a hybrid nanofluid and found that the capacity of the heat transfer rate of the hybrid nanofluid is greater than the regular nanofluid. The same state of heat transfer rate was obtained by Lund et al. [19,20] during the examination of a hybrid nanofluid. Yan et al. [21] also found similar results when they analyzed a hybrid nanofluid over the exponential surface with joule heating effects. Further, it is stated that “the skin friction coefficient, $f^{″}\left(0\right)$, enhances for the first solution when the suction S and magnetic M parameters are increased, while $f^{″}\left(0\right)$ reduces for the higher effect of the velocity slip factor, δ”. Some important articles about hybrid nanofluids for various effects can be accessed in these articles [22,23,24,25,26].

Mathematical analysis of the fluid flow problems for multiple solutions is important since these solutions cannot be seen experimentally [27]. In this regard, researchers claimed that these solutions exist because of the existence of non-linearity in equations of the fluid model and depend on the values of the applied parameters [28,29]. Many researchers have discussed the importance and applications of multiple solutions/branches. These solutions play an important role in developing the industry’s alternative flow option in an emergency. According to Khashi'ie et al. [30], “if the problem has non-unique solutions but the researchers manage to find one solution only, there is a probability that the solution is the lower branch solution (unstable/not real). This will lead to the misinterpretation of the flow and heat transfer characteristics”. According to Mishra and DebRoy [31], “multiple solutions have many important applications when these are related to heat transfer because the final qualities and structure of many products of material processing in the industries can be improved by the concept of multiple solutions”. Analysis of the stability is necessary to determine the stable branch when multiple branches occur. Many researchers have stated that only a stable branch has a physical significance which means that only a stable branch can be used in practical applications. Weidman [32] recently discussed the possibility that more than one solution could also be stable. According to him, “since the triple solutions appear to be upper branches (which cannot be continued to their lower branches), then those solutions will also be stable”. Therefore, multiple solutions/branches were considered in this study along with their stability analysis due to their important applications.

After evaluating the published literature, the motivation of this work is to examine the heat transfer characteristics of the rotating, steady, and 3D flow of the hybrid nanofluid. According to our best knowledge, no such study has been carried out for a hybrid nanofluid especially for the multiple solutions/branches.

2. Mathematical Description of the Problem

The flow of a three-dimensional hybrid nanofluid on an exponentially elastic shrinking surface is considered in a rotating frame of reference $x, y, z$. The velocity of the surface is$u_w=-a{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$l$}\right.}$, where $a$ is the characteristic velocity of the surface (See Figure 1). Momentum boundary layers of a hybrid nanofluid flow with energy equations without viscous dissipation and thermal radiation can be described as

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

(1)

$$u\frac{\partial u}{\partial x}+w\frac{\partial u}{\partial z}-2\tilde{\mathsf{\Omega }}v=\frac{{\mu }_{hnf}}{{\rho }_{hnf}}\frac{{\partial }^{2}u}{\partial z^2}$$

(2)

$$uu\frac{\partial v}{\partial x}+w\frac{\partial v}{\partial z}+2\tilde{\mathsf{\Omega }}u=\frac{{\mu }_{hnf}}{{\rho }_{hnf}}\frac{{\partial }^{2}u}{\partial z^2}$$

(3)

$$uu\frac{\partial T}{\partial x}+w\frac{\partial T}{\partial z}=\frac{k_{hnf}}{{\left(\rho c_p\right)}_{hnf}}\frac{{\partial }^{2}T}{\partial y^2}$$

(4)

The related boundary conditions (BCs) (2–5) are

$$\{\begin{array}{c}v=0 , u=u_w, w=w_0{e}^{\frac{x}{2l}}, T=T_w at z=0\\ u \to 0, v \to 0, T \to T_{\infty }, as z \to \infty \end{array}$$

(5)

where ${\left(\rho c_p\right)}_{hnf}{\mu }_{hnf},,k_{hnf}$ and$,{\rho }_{hnf}$ are the corresponding heat capacity, dynamic viscosity, thermal conductivity, and density of hybrid nanofluid. Moreover, subscript $hnf$ shows the thermophilic properties of hybrid nanofluid. Further, $w_0>0$ indicates the suction, and $w_0<0$ indicates the injection, and $\tilde{\mathsf{\Omega }}={\mathsf{\Omega }}_0{e}^{-\frac{x}{l}}$ is the local rotation parameter. The thermophysical properties are given in

Table 1. Thermophysical features of hybrid nanofluid.

Properties	Hybrid Nanofluid
Dynamic viscosity	${\mu }_{hnf}=\frac{{\mu }_f}{{\left(1-{\varphi }_{A{l}_2{O}_3}\right)}^{2.5}{\left(1-{\varphi }_{Cu}\right)}^{2.5}}$
Density	${\rho }_{hnf}=\left(1-{\varphi }_{Cu}\right)\left[\left(1-{\varphi }_{A{l}_2{O}_3}\right){\rho }_f+{\varphi }_{A{l}_2{O}_3}{\rho }_{A{l}_2{O}_3}\right]+{\varphi }_{Cu}{\rho }_{Cu}$
Thermal conductivity	$k_{hnf}=\frac{k_{Cu}+2k_{nf}-2{\varphi }_{Cu}\left(k_{nf}-k_{Cu}\right)}{k_{Cu}+2k_{nf}+{\varphi }_{Cu}\left(k_{nf}-k_{Cu}\right)}\times \left(k_{nf}\right)$ where $k_{nf}=\frac{k_{A{l}_2{O}_3}+2k_f-2{\varphi }_{A{l}_2{O}_3}\left(k_f-k_{A{l}_2{O}_3}\right)}{k_{A{l}_2{O}_3}+2k_f+{\varphi }_{A{l}_2{O}_3}\left(k_f-k_{A{l}_2{O}_3}\right)}\times \left(k_f\right)$
Heat capacity	${\left(\rho c_p\right)}_{hnf}=\left(1-{\varphi }_{Cu}\right)\left[\left(1-{\varphi }_{A{l}_2{O}_3}\right){\left(\rho c_p\right)}_f+{\varphi }_{A{l}_2{O}_3}{\left(\rho c_p\right)}_{A{l}_2{O}_3}\right]+{\varphi }_{Cu}{\left(\rho c_p\right)}_{Cu}$

and

Table 2. The thermophysical properties of nanofluid.

Fluids	P (kg/m3)	${\mathit{C}}_{\mathit{p}}(\mathrm{J}/\mathrm{kg}\mathrm{K})$	K (W/m K)
Alumina ($A{l}_2{O}_3$)	3970	765	40
Copper (Cu)	8933	385	400
Water ($H_{2}O$)	997.1	4179	0.613

.

We will employ similarity variable (6) in Equations (1)–(4) in order to obtain the similarity solutions

$$u=a{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$l$}\right.}{f}^{\prime }\left(\eta \right) ,v=a{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$l$}\right.}g\left(\eta \right), \theta \left(\eta \right)=\frac{\left(T-T_{\infty }\right)}{\left(T_w- T_{\infty }\right)}, \eta =z\sqrt{\frac{\vartheta a}{2l}} e^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$2l$}\right.}$$

(6)

Using the relationship between Equations (1) and (6), we obtain

$$w=\sqrt{\frac{\vartheta a}{2l}} e^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$2l$}\right.}\left\{f\left(\eta \right)+\eta f^{\prime }\left(\eta \right)\right\}$$

(7)

which leads to

$$f\left(0\right)=w_o\sqrt{\frac{2l}{\vartheta a}}=S$$

(8)

Substituting the stream function relationship with Equations (6)–(7) in Equations (2)–(5) yields

$$f^{‴}+{\xi }_1{\xi }_2\left[f{f}^{″}-2f^{{}^{\prime }2}+4\mathsf{\Omega }g\right]=0$$

(9)

$$g^{″}+{\xi }_1{\xi }_2\left[f{g}^{\prime }-2f^{\prime }g-4\mathsf{\Omega }{f}^{\prime }\right]=0$$

(10)

$$\frac{\left(k_{hnf}/k_f\right){\xi }_3}{Pr} {\theta }^{″}+f{\theta }^{\prime } =0$$

(11)

along with BCs

$$\{\begin{array}{c}f\left(0\right)=S, f^{\prime }\left(0\right)=-1, g\left(0\right)=0,\theta \left(0\right)=1\\ f^{\prime }\left(\eta \right) \to 0,g\left(\eta \right) \to 0 \theta \left(\eta \right)\to 0, as \eta \to \infty \end{array}$$

(12)

$$\{\begin{array}{c}{\xi }_1=\left\{\left(1-{\varphi }_{Cu}\right)\left[1-{\varphi }_{A{l}_2{O}_3}+{\varphi }_{A{l}_2{O}_3}\left(\raisebox{1ex}{${\rho }_{A{l}_2{O}_3}$}\!\left/ \!\raisebox{-1ex}{${\rho }_f$}\right.\right)\right]+{\varphi }_{Cu}\left(\raisebox{1ex}{${\rho }_{Cu}$}\!\left/ \!\raisebox{-1ex}{${\rho }_f$}\right.\right)\right\}\\ {\xi }_2={\left(1-{\varphi }_{Cu}\right)}^{2.5}{\left(1-{\varphi }_{A{l}_2{O}_3}\right)}^{2.5}\\ {\xi }_3=\frac{1}{\left\{\left(1-{\varphi }_{Cu}\right)\left[1-{\varphi }_{A{l}_2{O}_3}+{\varphi }_{A{l}_2{O}_3}\frac{{\left(\rho c_p\right)}_{A{l}_2{O}_3}}{{\left(\rho c_p\right)}_f}\right]+{\varphi }_{Cu}\frac{{\left(\rho c_p\right)}_{Cu}}{{\left(\rho c_p\right)}_f}\right\}}\end{array}$$

(13)

where prime represents the differentiation with respect to $\eta$, $Pr=\frac{\vartheta }{\alpha }$ is Prandtl, and $\mathsf{\Omega }=\frac{{\mathsf{\Omega }}_{0}l}{a}$ is the constant dimensionless rotation parameter.

Physical quantities are the skin friction coefficient and local Nusselt, which are expressed as

$$C_{fx}=\frac{{\mu }_{hnf}}{{\rho }_f{a}^2}\left(\frac{\partial u}{\partial z}\right)|z=0,C_{fy}=\frac{{\mu }_{hnf}}{{\rho }_f{w}_0^2}\left(\frac{\partial v}{\partial z}\right)|z=0,N{u}_x=-\frac{k_{hnf}}{k_f\left(T_w-T_{\infty }\right)}\left(\frac{\partial T}{\partial z}\right)|z=0$$

(14)

By substituting Equation (6) in Equation (14), the following is obtained

$$\begin{array}{l}2\sqrt{Re}{C}_{fx}=\frac{1}{{\left(1-{\varphi }_{A{l}_2{O}_3}\right)}^{2.5}{\left(1-{\varphi }_{Cu}\right)}^{2.5}}{f}^{″}\left(0\right);\\ 2\sqrt{Re}{C}_{fy}=\frac{1}{{\left(1-{\varphi }_{A{l}_2{O}_3}\right)}^{2.5}{\left(1-{\varphi }_{Cu}\right)}^{2.5}}{g}^{\prime }\left(0\right) \sqrt{\frac{Re}{2}}N{u}_x=-\frac{k_{hnf}}{k_f}{\theta }^{\prime }\left(0\right)\end{array}$$

(15)

where $R{e}_x=\frac{al}{\vartheta }{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$l$}\right.}$ is the local Reynold number.

3. Temporal Stability Analysis

The findings of the boundary layer problem (9–11) demonstrate that multiple branches occur. An analysis of stability is then carried out which has been performed by many researchers [33,34,35,36]. For the stability study, the unsteady forms of equations are supposed to be used. Henceforth, Equations (2)–(4) can be expressed as

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+w\frac{\partial u}{\partial z}-2\tilde{\mathsf{\Omega }}v=\frac{{\mu }_{hnf}}{{\rho }_{hnf}}\frac{{\partial }^{2}u}{\partial z^2}$$

(16)

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+w\frac{\partial v}{\partial z}+2\tilde{\mathsf{\Omega }}u=\frac{{\mu }_{hnf}}{{\rho }_{hnf}}\frac{{\partial }^{2}u}{\partial z^2}$$

(17)

$$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+w\frac{\partial T}{\partial z}=\frac{k_{hnf}}{{\left(\rho c_p\right)}_{hnf}}\frac{{\partial }^{2}T}{\partial y^2}$$

(18)

where $t$ indicates the time. By considering $t$ in terms of $\tau$, the current similarity variables (6–7) for the unsteady flow are as follows.

$$\{\begin{array}{c} u=a{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$l$}\right.}{f}_{\eta }\left(\eta ,\tau \right),v=a{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$l$}\right.}g\left(\eta ,\tau \right), \eta =z\sqrt{\frac{\vartheta a}{2l}} e^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$2l$}\right.} \\ w=\sqrt{\frac{\vartheta a}{2l}} e^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$2l$}\right.}\left\{f\left(\eta ,\tau \right)+\eta f_{\eta }\left(\eta ,\tau \right)\right\}, \tau =\frac{a}{2l}{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$l$}\right.}.t\\ \theta \left(\eta ,\tau \right)=\raisebox{1ex}{$\left(T-T_{\infty }\right)$}\!\left/ \!\raisebox{-1ex}{$\left(T_w-T_{\infty }\right)$}\right.\end{array}$$

(19)

By substituting Equation (19) into Equations (16)–(18), we obtain

$$f_{\eta \eta \eta }+{\xi }_1{\xi }_2\left[f{f}_{\eta \eta }-2{\left(f_{\eta }\right)}^2+4\mathsf{\Omega }g-f_{\tau \eta }\right]=0$$

(20)

$$g_{\eta \eta }+{\xi }_1{\xi }_2\left[f{g}_{\eta }-2f_{\eta }g-4\mathsf{\Omega }{f}_{\eta }-g_{\tau }\right]=0$$

(21)

$$\frac{\left(k_{hnf}/k_f\right){\xi }_3}{Pr} {\theta }_{\eta \eta }+{\theta }_{\eta }f-{\theta }_{\tau }=0$$

(22)

While BCs (12) can be:

$$\{\begin{array}{c}f\left(0,\tau \right)=S, f_{\eta }\left(0,\tau \right)=-1, g\left(0,\tau \right)=0,\theta \left(0,\tau \right)=1\\ f_{\eta }\left(\eta ,\tau \right) \to 0,g\left(\eta ,\tau \right) \to 0 \theta \left(\eta ,\tau \right)\to 0, as \eta \to \infty \end{array}$$

(23)

According to Weidman et al. [37], “the stability of the steady flow solutions $f\left(\eta \right)=f_0\left(\eta \right),g\left(\eta \right)=g_0\left(\eta \right)$ and $\theta \left(\eta \right)={\theta }_0\left(\eta \right)$ are identified by writing $F\left(\eta ,\tau \right), G\left(\eta ,\tau \right)$, and $H\left(\eta ,\tau \right)$” as follows

$$\begin{gathered} f\left(\eta ,\tau \right)=f_0\left(\eta \right)+e^{-\epsilon \tau }F\left(\eta ,\tau \right), g\left(\eta ,\tau \right)=g_0\left(\eta \right)+e^{-\epsilon \tau }G\left(\eta ,\tau \right), \theta \left(\eta ,\tau \right) \\ ={\theta }_0\left(\eta \right)+e^{-\epsilon \tau }H\left(\eta ,\tau \right) \end{gathered}$$

(24)

where $\epsilon$ is the unknown eigenvalue parameter, while functions $F\left(\eta ,\tau \right), G\left(\eta ,\tau \right)$, Hη, and $\left(\tau \right)$ are small relative to$f\left(\eta \right)=f_0\left(\eta \right),g\left(\eta \right)=g_0\left(\eta \right)$, and $\theta \left(\eta \right)={\theta }_0\left(\eta \right)$. The following system of equations is obtained by putting Equation (24) into Equations (20)–(23).

$$F_{\eta \eta \eta }+{\xi }_1{\xi }_2\left[f_0{F}_{\eta \eta }+F{\left(f_0\right)}_{\eta \eta }-4{\left(f_0\right)}_{\eta }{F}_{\eta }+4\mathsf{\Omega }G+\epsilon F_{\eta }\right]=0$$

(25)

$$G_{\eta \eta }+{\xi }_1{\xi }_2\left[f_0{G}_{\eta }+F{\left(g_0\right)}_{\eta }-2{\left(f_0\right)}_{\eta }G-2F_{\eta }{g}_0-4\mathsf{\Omega }{F}_{\eta }+\epsilon G\right]=0$$

(26)

$$\frac{\left(k_{hnf}/k_f\right){\xi }_3}{Pr}{H}_{\eta \eta }+f_0{H}_{\eta }+F{\left({\theta }_0\right)}_{\eta }+\epsilon H=0$$

(27)

with the following BCs

$$\{\begin{array}{c}F\left(0,\tau \right)=0, F_{\eta }\left(0,\tau \right)=0, G\left(0,\tau \right)=0,H\left(0,\tau \right)=0\\ F_{\eta }\left(\eta ,\tau \right) \to 0,G\left(\eta ,\tau \right) \to 0, H\left(\eta ,\tau \right)\to 0, as \eta \to \infty \end{array}$$

(28)

The stabilization of heat transfer solutions and steady-state flow solutions $f_0\left(\eta \right), g_0\left(\eta \right)$, and ${\theta }_0\left(\eta \right)$ can be obtained by setting $\tau \to 0.$ Therefore, the functions $F\left(\eta ,\tau \right)=F_0\left(\eta \right)$, $G\left(\eta ,\tau \right)=G_0\left(\eta \right)$, and $H\left(\eta ,\tau \right)=H_0\left(\eta \right)$ can be written in Equations (25)–(27). Thus, the following problem of linearized eigenvalue can be expressed as:

$$F_0^{‴}+{\xi }_1{\xi }_2\left\{f_0{F}_0^{″}+F_0{f}_0^{″}-4f_0^{\prime }{F}_0^{\prime }+4\mathsf{\Omega }{G}_0+\epsilon F_0^{\prime }\right\}=0$$

(29)

$$G_0^{″}+{\xi }_1{\xi }_2\left[G_0^{\prime }{f}_0+g_0^{\prime }{F}_0-2\left(f_0^{\prime }{G}_0+F_0^{\prime }{g}_0\right)-4\mathsf{\Omega }{F}_0^{\prime }+\epsilon G_0\right]=0$$

(30)

$$\frac{\left(k_{hnf}/k_f\right){\xi }_3}{Pr}{H}_0^{″}+{\theta }_0^{\prime }{F}_0+H_0^{\prime }{f}_0+\epsilon H_0=0$$

(31)

subject to the following BCs:

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

(32)

Note that the smallest eigenvalue ($\epsilon$) can be determined by easing the boundary condition [38,39]. In this analysis, the condition $F_0^{\prime }\left(\eta \right)\to 0$ was relaxed and Equations (29)–(31) were solved along with a new relaxed BC $F_0^{″}\left(0\right)=1$ for a fixed value of the applied parameters.

4. Numerical Method

The three-stage Labatto III-A method is adopted to solve the above system of Equations (9)–(12) numerically with the help of a bvp4c solver in the MATLAB software. Shampine et al. [40] provides a detailed explanation of this method and how it works in MATLAB. We kept the tolerance at ${10}^{-5}$to obtain good accuracy in the solutions. According to Raza et al. [41] and Lund et al. [42], “this collocation polynomial and formula offers a $C^1$ continuous solution in which mesh error control and selection are created on the residual of the continuous solution. The tolerance of relative error is fixed ${10}^{-5}$ for the current problem. The suitable mesh determination is created and returned in the field sol.x. The bvp4c returns solution, called as sol.y., as a construction. In any case, values of the solution are gotten from the array named sol.y relating to the field sol.x”. To determine the stable solution, the three-stage Labatto III-A method is also used to obtain the values of the smallest eigenvalue. For a better understanding, the algorithm of the method is illustrated in Figure 2.

5. Results and Discussion

The system of nonlinear ODEs (9–11) along with BCs (12) was solved numerically by employing the three-stage Labatto III-A method. We have kept fixed $Pr=6.2$ for the water-based hybrid nanofluid at room temperature 25 °C. The Using Keller box method, the obtained results were compared with the published results of Rosali et al. [43] as shown in

Table 3. Comparison of $f^{″}\left(0\right)$ and $g^{\prime }\left(0\right)$ for various values of $\mathsf{\Omega }$ when $=2, Pr=1, {\varphi }_{Cu}={\varphi }_{A{l}_2{O}_3}=0$.

	Results of [43]		Present Results
$\mathsf{\Omega }$	$f^{″}\left(0\right)$	$g^{\prime }\left(0\right)$	$f^{″}\left(0\right)$	$g^{\prime }\left(0\right)$	$-{\theta }^{\prime }\left(0\right)$
0.2	1.1449077	0.76030466	1.14491145	0.76028572	1.65895348
0.5	1.7030846	1.06450099	1.70307504	1.06449792	1.72967817
1	2.2015332	1.45971480	2.20152305	1.45971069	1.77530624
2	2.8466686	2.03340224	2.84665295	2.03339910	1.81789928
5	4.0629311	3.18470073	4.06290502	3.18469911	1.86820205
10	5.4008517	4.48862880	5.40081392	4.48862871	1.90024712

for the coefficients of the skin friction $f^{″}\left(0\right)$ and $g^{\prime }\left(0\right)$ for numerous values of the local rotating parameter $\mathsf{\Omega }$ and found to be in excellent agreement. Therefore, it can be concluded that the current numerical technique can be employed with considerable confidence to solve the considered problem.

Table 4. Values of $f^{″}\left(0\right)$, $g^{\prime }\left(0\right)$ and $-{\theta }^{\prime }\left(0\right)$ for the various values of ${\varphi }_{Cu},{\varphi }_{A{l}_2{O}_3}$ when $Pr=6.2,S=2.5,\mathsf{\Omega }=0.01$.

			Upper			Lower
${\varphi }_{A{l}_2{O}_3}$	${\varphi }_{Cu}$	$f^{″}\left(0\right)$	$g^{\prime }\left(0\right)$	$-{\theta }^{\prime }\left(0\right)$	$f^{″}\left(0\right)$	$g^{\prime }\left(0\right)$	$-{\theta }^{\prime }\left(0\right)$
0	0.01	1.803107	0.026679	14.648604	−0.018593	0.250665	14.599120
	0.05	2.304795	0.023603	12.924215	−0.498357	0.454046	12.840883
	0.1	2.710236	0.022151	11.100902	−0.977529	0.682433	10.978038
0.01	0.01	1.808890	0.026631	14.188081	−0.023519	0.252548	14.136747
	0.05	2.287739	0.023678	2.5184567	−0.480112	0.445739	12.433579
	0.1	2.672587	0.022263	10.752889	−0.929227	0.658664	10.629077
0.1	0.01	1.735351	0.027287	10.736454	−0.0382667	0.229348	10.673921
	0.05	2.045129	0.024944	9.4760101	−0.2360146	0.338429	9.3837847
	0.1	2.275158	0.023735	8.1418098	−0.4667501	0.439678	8.0172755

was constructed for the values of the $f^{″}\left(0\right)$ $g^{\prime }\left(0\right)$ and $-{\theta }^{\prime }\left(0\right)$ of hybrid nanofluid.

The effects of ${\varphi }_{Cu}$ on the coefficients of skin friction $f^{″}\left(0\right)$, $g^{\prime }\left(0\right)$, and the heat transfer rate $-{\theta }^{\prime }\left(0\right)$ against various values of the suction parameter $S$ are given in Figure 3, Figure 4 and Figure 5. Here, we focus solely on multiple branches. The non-uniqueness of the branches is only possible when $\mathsf{\Omega }=0.01$ (refer to Figure 6). Furthermore, dual branches for the Equations (9)–(11) exist when $S\le S_c$ where $S_c$ is the critical value. No solution exists when $S>S_c$. It is observed that the corresponding critical values of ${\varphi }_{Cu}=0.01, 0.05, 0.1$ are $S_c=1.3914, 1.3249, 1.2793$, respectively.

The skin friction increases (decreases) by increasing the copper volume fraction in the upper (lower) branch. It is also examined that when the volume fraction of the copper enhances, the separation of the boundary layer expands. The heat transfer rate reduces when ${\varphi }_{Cu}$ increases in both branches, while it is increasing the function of the suction parameter.

Figure 6, Figure 7 and Figure 8 were plotted to demonstrate the effect of $\mathsf{\Omega }$ on $f^{″}\left(0\right)$,$g^{\prime }\left(0\right)$, and $-{\theta }^{\prime }\left(0\right)$ against the fixed values of $S$. It was observed that skin friction coefficient $f^{″}\left(0\right)$ and heat transfer $-{\theta }^{\prime }\left(0\right)$ are increasing functions of the rotating ($\mathsf{\Omega }$) parameter when suction $S$ increases. It is also shown that $f^{″}\left(0\right)$ increases for the higher values of rotational parameter $\mathsf{\Omega }$ in both positive and negative sides. Moreover, in the case of $\mathsf{\Omega }>0$, $g^{\prime }\left(0\right)$ increases when $\mathsf{\Omega }$ increases. Furthermore, it is revealed from Figure 7 that the skin friction increases with the decelerated flow. (i.e, $\mathsf{\Omega }<0$) and decreases with the accelerated flow (i.e, $\mathsf{\Omega }>0$). On the other hand, no solution is found when $\mathsf{\Omega }=0$ for the fixed value of $S=1$. Furthermore, the symmetrical behavior of the branches is shown in these figures.

Figure 9 was plotted to examine the effects of $\mathsf{\Omega }$ on the hybrid nanofluid velocity $f^{\prime }\left(\eta \right)$. It was detected that the velocity of the hybrid nanofluid declines as the rotation ($\mathsf{\Omega }$) parameter is increased. Moreover, no oscillation behavior is found in $f^{\prime }\left(\eta \right)$ for the higher values of $\mathsf{\Omega }$. It happened due to various effects, such as the effects of the shrinking parameter, volume fraction, and suction.

Figure 10 was plotted to examine the effects of the rotation parameter $\mathsf{\Omega }$ on the hybrid nanofluid velocity $g\left(\eta \right)$. When $\mathsf{\Omega }$ is increased, the velocity of hybrid nanofluid contains duality in the behavior. For the negative and positive values of the rotation ($\mathsf{\Omega })$ parameter, the behavior of the velocity profile was found to have the same behavior. Physically, this indicates that there is a symmetrical solution to the hybrid nanofluid problem.

The values of the smallest eigenvalues $\epsilon$ for different values of $S$ and ${\varphi }_{Cu}$ are shown in

Table 5. Smallest eigenvalue $\epsilon$ for different values of S and ${\varphi }_{Cu}$ when $\mathsf{\Omega }=0.01; Pr=6.2; {\varphi }_{A{l}_2{O}_3}=0.1$.

${\mathit{\varphi }}_{\mathit{C}\mathit{u}}$	$\mathit{S}$	$\mathit{\epsilon }$
		Upper branch	Lower branch
0.01	1.4	0.0527	−0.02652
	1.6	0.26126	−0.13272
	1.8	0.56249	−0.63379
0.05	1.5	0.15122	−0.17922
	1.7	0.42821	−0.50293
	1.9	0.79456	−0.73434
0.1	1.5	0.26055	−0.27788
	1.7	0.49440	−0.53679
	1.9	0.84495	−0.80766

. The positive eigenvalue causes the initial decay of disturbance and thus stabilizes the flow. In contrast, the negative results of the smallest eigenvalue show that the flow is unstable. Table 5 shows that $\epsilon$ is positive for the upper branch, whereas $\epsilon$ is negative for the lower branch.

6. Conclusions

In this study, the flow of rotating, steady, and three-dimensional heat transfer of a hybrid nanofluid on a penetrable exponential shrinking surface together with the suction effect were investigated numerically. The governing PDEs have been converted to a system of ODEs using the suitable exponential similarity transformation. The three-stage Labatto III-A technique was then implemented for the solving of the system of ODEs. Numerical results indicate that the current outcomes of $f^{″}\left(0\right)$ and $g^{\prime }\left(0\right)$ are in good agreement with the results previously published. The point-wise conclusions are the following:

• In comparison to a viscous fluid, the heat transfer rate of the hybrid nanofluid is better in the attendance of hybrid nanoparticles.
• Two branches exist in the specific ranges of physical parameters.
• The upper branch remains stable while the lower branch is unstable.
• Rate of heat transfer upsurges for the advanced values of the mass suction in both branches.
• Hybrid nanofluid has a symmetrical solution.