1. Introduction
The thermal performance of fluids employed for heat transport analysis is of main concern in various thermal processes used in industry. In recent years, several researchers put their efforts into improving the thermal performance of these processes by employing different techniques. Correspondingly, it was shown that the thermal features of ordinary fluids were substantially changed with the addition of nanosized particles of single type, known as nanoparticles. Such nanoparticles could be made up of oxides, carbides, carbon nanotubes, and metals. Nanofluids are a stable mixture of base fluid and suspended nanoparticles (up to 100 nm). Engine oil, water, and ethylene glycol may be taken as base fluids. As expected, nanofluids displayed higher heat dissipation performance. In addition, they have shown greater increases in their thermal physical characteristics, like thermal conductivity, etc. Treatment of cancer, medicine, transformer cooling, solar collectors, heat exchangers, nuclear reactor cooling, freezers, electronics cooling, and automobiles are just a few examples of nanofluid applications. Due to the successful performance of nanofluids in the heat transport mechanism, numerous studies have reported on the heat transfer of nanofluids. Choi [1] was the first one to spotlight the notion and implications of using nanoparticles in a conventional fluid to improve heat transfer. After that, Buongiorno [2] and Tiwari and Das [3] proposed two common hydrodynamic models for nanofluids. Buongiorno's model considered the impacts of Brownian motion and thermophoresis. On the other hand, the Tiwari and Das model explored the nature of nanofluids by volume fraction considering the nanoparticles. Several authors, including Kuznetsov and Nield [4], Makinde and Aziz [5], Mustafa et al [6], Khan and Aziz [7], and Hashim and Khan [8], investigated the various physical characteristics of nanofluids using Buongiorno's model. On the other hand, the single-phase model proposed by Tiwari and Das [3] has been applied by various authors. Rashad et al [9] utilized the numerical technique to explore the mixed convection flow of copper–water nanofluid in a rectangular cavity within a porous medium. The natural convection flow of an Al2O3/water nanofluid along with heat transport in an L-shaped enclosure has been presented by Mohebbi and Rashidi [10]. Later, Bhatti et al [11] numerically investigated the 3D unsteady flow of a nanofluid with gyrotactic microorganisms driven by cylindrical geometry. Many researchers have explored the flow and heat transport properties in various geometries using nanofluids in recent years due to the effectiveness of nanofluids, including Hashim et al [12], Hamid et al [13], Hafeez et al [14], etc.
These days, the research community focuses on radiative heat transfer and flow processes in energy conversion systems that run at high temperatures, owing to their remarkable performance in a variety of fields of science, including satellites, missiles, various aircraft propulsion devices, gas turbines, and nuclear power plants. For a device operating at an above-average temperature where radiation from heated walls and working fluid is different, the effect of thermal radiation is highly important in the transfer of heat flow. Heat transfer features under the influence of thermal radiation have been investigated by several researchers. The influences of buoyancy force and thermal radiation on stagnation point flow past a stretching sheet were examined by Pal [15]. Furthermore, Bidin and Nazar [16] numerically investigated the effects of thermal radiation on flow and heat transport analysis due to an exponentially extending surface. Later, Aziz [17] conducted a numerical study on the flow and heat transport mechanism for viscous fluid flow over an unsteady stretching surface. Pal and Mondal [18] analysed non-Darcy flow over a stretching plate incorporating thermal radiation effects. Dogonchi and Ganji [19] addressed the heat transport characteristics of a nanofluid past a stretching surface in the presence of thermal radiation. Lin et al [20] explored the flow and heat transport characteristics of copper–water nanofluid subject to thermal radiation and nanoparticle shape factor effects. Waqas et al [21] presented a numerical study to address the flow of non-Newtonian Carreau fluid subject to thermal radiation using revised nanofluid model. Later, Sheikholeslami et al [22] employed the control volume-based finite element numerical method to study the flow of a nanofluid through a wavy chamber with thermal radiation impacts. The closed-form solutions for two-dimensional flow and heat transfer analysis in the presence of thermal radiation past a vertical plate have been examined by Turkyilmazoglu [23].
Stratification occurs because of temperature gradients or the mixture of multiple fluids having different densities in engineering and industrial mechanisms. Exploring the mechanism of thermal stratification during flow and heat transport in nanofluids is of tremendous physical interest. This process is essential in the disciplines of lake thermo-hydraulics, salinity and thermal stratification mechanisms in oceans, heat rejection into the surrounding environment via rivers, agriculture fields, volcanic flows, and industries such as reservoirs, industrial food and salinity, atmosphere involving heterogeneous mixtures, and groundwater reservoirs. As a pioneer, Yang et al [24] investigated the free convective flow of a thermally stratified fluid due to a non-isothermal plate. Ishak et al [25] examined the time-independent mixed convection flow through a stable stratified medium near a vertical flat plate. Cheng [26] employed the cubic spline collocation method to explore the impacts of double stratification on natural convection flow of a non-Newtonian fluid near a vertical wavy surface. Mukhopadhyay and Ishak [27] presented a numerical study to discuss the mixed convection axisymmetric flow of thermally stratified viscous fluid over a stretching cylinder. Mishra et al [28] conducted a numerical study to probe a steady flow due to a vertical surface subject to double stratified micropolar fluid. An identical study of nanofluid flow and heat transport with thermal stratification has been presented by Abbasi et al [29], Hayat et al [30], Eswaramoorthi et al [31], Jabeen et al [32].
According to the literature review, there have been several investigations on the flow and transfer of hybrid nanofluids using various mechanisms. In most of these works, the authors have studied the hybrid nanofluid heat transport features by considering constant wall and free stream temperatures and computed single solutions by using analytical and numerical techniques. However, in a variety of real-world circumstances, these temperatures do not remain constant, and we must treat them as a function of space and time variables. Moreover, the energy transport phenomenon of a hybrid nanofluid driven by a stretching/shrinking geometry with variable temperatures has various realistic industrial, engineering and biomechanical applications, like polymer technology, blood flow, treatment of several diseases, metallurgical processes, and annealing and thinning of wires. As a result of these applications in numerous disciplines of science and technology, it is worthwhile to discuss and explore thermally stratified flow of hybrid nanofluids with thermal radiation.
As per the authors' knowledge and based on the open literature review, it is noticed that multiple numerical solutions for thermally stratified flow of hybrid nanofluid driven by a shrinking surface have not been reported yet. The core novelty of the current study is to perform a numerical simulation to predict the multiple solutions for thermally stratified flow of ${\rm{Cu}}-{{\rm{Al}}}_{2}{{\rm{O}}}_{3}-$ hybrid nanomaterials along with heat transport analysis in the presence of thermal radiation and slip mechanism. In this research, the authors formulated the problem of two-dimensional time-dependent magnetohydrodynamic (MHD) hybrid nanofluid flow over a flat sheet with the help of conservation laws and Boussinesq approximations in the form of partial differential equations. Moreover, the current investigation employs an efficient and versatile numerical method, the bvp4c routine in MATLAB, to acquire the problem solutions. The outcomes are compared to those obtained without hybrid nanofluids as well as those obtained with nanofluids.
2. Mathematical model
2.1. Formulation and basic equations
As shown in figure 1, the considered physical situation involves unsteady, two-dimensional, and incompressible flow of hybrid nanofluids past a flat surface with stretching and shrinking characteristics. During the mathematical formations, the following assumptions have been kept in mind:
| 1. The Cartesian coordinate system $(x,\,y)\,$ is chosen as a frame of reference. | |
| 2. The single-phase model is utilized to describe the hybrid nanofluids while the Rosseland approximation model is taken for the radiative heat transport. | |
| 3. The hybrid nanofluids are composed of copper $\left({\rm{Cu}}\right)\,$ and aluminium oxide $\left({{\rm{Al}}}_{{\rm{2}}}{{\rm{O}}}_{{\rm{3}}}\right)\,$ nanoparticles with water as the base fluid. | |
| 4. The flat surface is being stretched or shrunk with velocity ${u}_{w}=\tfrac{ax}{1-\alpha t},$ where $a$ signifies the shrinking/stretching rate on the $x-$ axis. | |
| 5. The impacts of a variable magnetic field $B(t)=\tfrac{{B}_{o}}{\sqrt{1-\alpha t}}\,$ are taken in the vertical direction to the flow by neglecting the induced magnetic field. | |
| 6. Thermal stratification impacts are imposed on the flow fields for which the wall and ambient temperatures are taken to be variable, with the values ${T}_{w}={T}_{0}+\tfrac{{m}_{1}x}{1-\alpha t}$ and ${T}_{\infty }={T}_{0}+\tfrac{{m}_{2}x}{1-\alpha t}$, respectively. |
Figure 1. Flow geometry and coordinate system. |
Therefore, in view of the Boussinesq approximations and above stated restrictions, the basic conservation equations for MHD unsteady flow of hybrid nanofluids subject to the Tiwari and Das [3] model can be acknowledged as (see Fang et al [33], Rohini et al [34], Devi and Devi [35], Ismail et al [36]):
Continuity equation:
$\begin{eqnarray}\displaystyle \frac{\partial u}{\partial x}+\displaystyle \frac{\partial v}{\partial y}=0,\,\end{eqnarray}$
Momentum equation: $\begin{eqnarray}\begin{array}{l}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial \bar{y}}=\frac{\partial {u}_{\infty }}{\partial t}+{u}_{\infty }\frac{\partial {u}_{\infty }}{\partial x}\\ \,+\,\frac{{\mu }_{\mathrm{hnf}}}{{\rho }_{\mathrm{hnf}}}\frac{{\partial }^{2}u}{\partial {y}^{2}}+\frac{{\sigma }_{f}{B}^{2}(t)}{{\rho }_{\mathrm{hnf}}}\left({u}_{\infty }-u\right),\,\end{array}\end{eqnarray}$
Energy equation: $\begin{eqnarray}\begin{array}{l}\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial \bar{y}}=\frac{{k}_{\mathrm{hnf}}}{{\left(\rho {C}_{p}\right)}_{\mathrm{hnf}}}\frac{{\partial }^{2}T}{\partial {y}^{2}}\\ \,-\,\frac{1}{{\left(\rho {C}_{p}\right)}_{\mathrm{hnf}}}\frac{\partial {q}_{r}}{\partial y}+\frac{{Q}^{\ast }}{{\left(\rho {C}_{p}\right)}_{\mathrm{hnf}}}\left(T-{T}_{\infty }\right),\end{array}\end{eqnarray}$
where $\left(u,v\right)$ represents the velocity components along the $(x,\,y)\,-$axis, $\left({\rho }_{\mathrm{hnf}},{\mu }_{\mathrm{hnf}},{\left({c}_{p}\right)}_{\mathrm{hnf}},{k}_{\mathrm{hnf}}\right)$ means the density, viscosity, specific heat and thermal conductivity of hybrid nanofluids, ${u}_{\infty }$ signifies the free stream velocity, and ${B}_{0}\,$ is the magnitude of the applied magnetic field. Furthermore, the radiative heat flux is denoted by ${q}_{r,\mathrm{rad}}$ and ${Q}^{\ast }=\tfrac{{Q}_{0}}{1-\alpha t},$ means the heat generation/absorption coefficient.The following boundary conditions are put on the surface walls as well as on the free stream:4 ) we have imposed the velocity-slip boundary condition along with suction or injection at the surface wall. Physically, the velocity-slip mechanism occurs due to a difference between the wall velocity and the velocity of the fluid particles adjacent to the wall. Mathematically, the velocity-slip is taken to be proportional to the local shear stress. Slip flow is typically specified by a nonphysical quantity called the slip length ${L}_{1}^{\ast }.$ The second condition represents the suction or injection phenomenon at the wall which is described by the mass transfer velocity ${v}_{w}\left(x,t\right).$ For temperature fields, we have considered the variables temperatures ${T}_{w}$ and ${T}_{\infty }$ along the wall and away from the wall. The physical properties of the hybrid nanofluids and the base fluid are already provided in table 1. In addition, the theoretical models for calculating these physical properties of hybrid nanofluids are presented by the following expressions:
$\begin{eqnarray}{\rm{At}}\,y=0;\,u={u}_{w}+{L}_{1}^{\ast }{\nu }_{f}\displaystyle \frac{\partial u}{\partial y},\,v={v}_{w}\left(x,t\right),\,T={T}_{w},\,\end{eqnarray}$
$\begin{eqnarray}{\rm{As}}\,y\to \infty ;\,u\to {u}_{\infty }\left(=\displaystyle \frac{bx}{1-\alpha t}\right),\,T\to {T}_{\infty },\end{eqnarray}$
where ${L}_{1}^{\ast }={L}_{1}\sqrt{1-\alpha t}\,$ represents the velocity-slip factor. Here, in equation ( $\begin{eqnarray*}{\mu }_{\mathrm{hnf}}=\frac{{\mu }_{f}}{{\left(1-{\phi }_{1}\right)}^{2.5}{\left(1-{\phi }_{2}\right)}^{2.5}},\,\end{eqnarray*}$
$\begin{eqnarray*}{\rho }_{\mathrm{hnf}}=\left(1-{\phi }_{2}\right)\,\left[\left(1-{\phi }_{1}\right)\,{\rho }_{f}+{\phi }_{1}{\rho }_{s1}\right]+{\phi }_{2}{\rho }_{s2},\,\end{eqnarray*}$
$\begin{eqnarray*}{\left(\rho {C}_{p}\right)}_{\mathrm{hnf}}=\left(1-{\phi }_{2}\right)\,\left[\left(1-{\phi }_{1}\right)\,{\left(\rho {C}_{p}\right)}_{f}+{\phi }_{1}{\left(\rho {C}_{p}\right)}_{s1}\right]+{\phi }_{2}{\left(\rho {C}_{p}\right)}_{s2},\,\end{eqnarray*}$
$\begin{eqnarray*}{k}_{\mathrm{hnf}}=\frac{{k}_{s2}+2{k}_{\mathrm{nf}}-2{\phi }_{2}\left({k}_{\mathrm{nf}}-{k}_{s2}\right)}{{k}_{s2}+2{k}_{\mathrm{nf}}+{\phi }_{2}\left({k}_{\mathrm{nf}}-{k}_{s2}\right)}{k}_{\mathrm{nf}},\end{eqnarray*}$
$\begin{eqnarray}{k}_{\mathrm{nf}}=\frac{{k}_{s1}+2{k}_{f}-2{\phi }_{1}\left({k}_{f}-{k}_{s1}\right)}{{k}_{s2}+2k+{\phi }_{1}\left({k}_{f}-{k}_{s1}\right)}{k}_{f}.\end{eqnarray}$
In the above equation, the subscript $ {``} \mathrm{hnf} {''} $ represents the hybrid nanofluids, while $ {``} f {''} $ denotes the base fluid.Table 1. Thermo-physical characteristics of hybrid nanofluids and base fluid (see Khanafer et al [37], Oztop and Abu-Nada [38]). |
| Physical properties | Base fluid | Nanoparticles | |
|---|---|---|---|
| ${{\rm{H}}}_{{\rm{2}}}{\rm{O}}$ | ${\rm{Cu}}$ | ${{\rm{Al}}}_{{\rm{2}}}{{\rm{O}}}_{{\rm{3}}}$ | |
| $\rho \,\left(\mathrm{kg}/{{\rm{m}}}^{3}\right)$ | 997.1 | 8933 | 3970 |
| ${c}_{p}\,\left({\rm{J}}/(\mathrm{kg}\,{\rm{K}})\right)$ | 4179 | 385 | 765 |
| $k\,\left({\rm{W}}/({\rm{m}}\,{\rm{K}})\right)$ | 0.613 | 400 | 40 |
The radiative heat flux ${q}_{r}\,$ in equation (3 ) is depicted as:7 ) into equation (3 ) , we get
$\begin{eqnarray}{q}_{r}=-\displaystyle \frac{4{\sigma }_{1}^{\ast }}{3{k}^{\ast }}\displaystyle \frac{\partial {T}^{4}}{\partial y}=-\displaystyle \frac{16{\sigma }^{\ast }}{3{k}^{\ast }}{T}^{3}\displaystyle \frac{\partial T}{\partial y},\,\end{eqnarray}$
where $\left({\sigma }_{1}^{\ast },{k}^{* }\right)$ signifies the Stefan–Boltzmann constant and coefficient of mean absorption. Now, plugging ${q}_{r}\,$ from equation ( $\begin{eqnarray}\begin{array}{l}\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\frac{{k}_{\mathrm{hnf}}}{{\left(\rho {C}_{p}\right)}_{\mathrm{hnf}}}\frac{{\partial }^{2}T}{\partial {y}^{2}}\\ \,+\,\frac{16{\sigma }^{\ast }{T}_{\infty }^{3}}{{\left(\rho {C}_{p}\right)}_{\mathrm{hnf}}3{k}^{\ast }}\frac{{\partial }^{2}T}{\partial {y}^{2}}+\frac{{Q}^{\ast }}{{\left(\rho {C}_{p}\right)}_{\mathrm{hnf}}}\left(T-{T}_{\infty }\right).\end{array}\end{eqnarray}$
The dimensionless form of the modelled problem is obtained by utilizing the following subsequent transformations: $\begin{eqnarray}\begin{array}{l}\psi =x\sqrt{\displaystyle \frac{b{\nu }_{f}}{1-\alpha t}}f\left(\eta \right),\,\,\theta \left(\eta \right)=\displaystyle \frac{T-{T}_{\infty }}{{T}_{w}-{T}_{o}},\,\\ \eta =y\sqrt{\displaystyle \frac{b}{{\nu }_{f}(1-\alpha t)}}.\,\end{array}\end{eqnarray}$
Substituting equation (9 ) into equations (2 ) and (8 ), the following converted system of ordinary differential equations is derived:
$\begin{eqnarray}\begin{array}{l}{f}^{\prime\prime\prime} +{\left(1-{\phi }_{1}\right)}^{2.5}{\left(1-{\phi }_{2}\right)}^{2.5}\\ \,\times \left\{\left(1-{\phi }_{2}\right)\,\left[1-{\phi }_{1}+{\phi }_{1}\left(\displaystyle \frac{{\rho }_{s1}}{{\rho }_{f}}\right)\right]+{\phi }_{2}\left(\displaystyle \frac{{\rho }_{s2}}{{\rho }_{f}}\right)\right\}\\ \,\times \left[ff^{\prime\prime} +1-{f^{\prime} }^{2}-\beta \left(f^{\prime} +\displaystyle \frac{\eta }{2}f^{\prime\prime} -1\right)\right]\\ \,+M{\left(1-{\phi }_{1}\right)}^{2.5}{\left(1-{\phi }_{2}\right)}^{2.5}\left(1-f^{\prime} \right)=0,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\frac{1}{\Pr }\left(\frac{{k}_{\mathrm{hnf}}}{{k}_{f}}+\frac{4}{3}R\right)\,\theta ^{\prime\prime} +\left\{\left(1-{\phi }_{2}\right)\,\left[1-{\phi }_{1}+{\phi }_{1}\frac{{\left(\rho {C}_{p}\right)}_{s1}}{{\left(\rho {C}_{p}\right)}_{f}}\right]\right.\\ \left.\,+\,{\phi }_{2}\frac{{\left(\rho {C}_{p}\right)}_{s2}}{{\left(\rho {C}_{p}\right)}_{f}}\right\}\,\left[\begin{array}{l}f\theta ^{\prime} -f^{\prime} \left(\theta +{s}_{t}\right)\\ -\beta \left(\theta +\frac{\eta }{2}\theta ^{\prime} +{s}_{t}\right)\end{array}\right]+\delta \theta =0,\end{array}\end{eqnarray}$
subject to the conditions $\begin{eqnarray}\begin{array}{l}f(0)=S,f^{'}(0)=\lambda +{\alpha }_{1}f^{''} (0),\,\theta (0)=1-{s}_{t},\,\,\\ f^{'}(\eta )\to 1,\theta \left(\eta \right)\to 0\,\mathrm{as}\,\eta \to \infty .\end{array}\end{eqnarray}$
The wall mass transfer velocity (or fluid suction velocity) becomes ${v}_{w}\left(0,t\right)=-\sqrt{\tfrac{b{\nu }_{f}}{1-\alpha t}}f\left(0\right)=S,$ where $S$ is a constant that specifies the wall mass transfer parameter, with $S\gt 0$ indicating suction, $S\lt 0$ indicating injection and $S=0$ indicating impermeability.The other dimensionless parameters are written as
$\begin{eqnarray}\begin{array}{l}{\alpha }_{1}={L}_{1}\sqrt{b{\nu }_{f}},M=\frac{{\sigma }_{f}{B}_{o}^{2}}{b{\rho }_{f}},\lambda =\frac{a}{b},\beta =\frac{\alpha }{b},\\ Pr =\frac{{v}_{f}}{{\alpha }_{f}},\delta =\frac{Q}{b{\left(\rho {C}_{p}\right)}_{f}},{s}_{t}=\frac{{m}_{2}}{{m}_{1}},R=\frac{4\sigma * {T}_{\infty }^{3}}{{k}_{f}k* }.\end{array}\end{eqnarray}$
In this analysis, the variables of physical curiosity are the drag force and the heat transfer rate, which are written in their dimensionless form as follows:
$\begin{eqnarray}\begin{array}{l}{C}_{f}{{Re}}_{x}^{1/2}=\frac{1}{{\left(1-{\phi }_{1}\right)}^{2.5}{\left(1-{\phi }_{2}\right)}^{2.5}}{f}^{{\prime\prime} }\left(0\right),N{u}_{x}{{{Re}}_{x}}^{-1/2}\\ \,\,\,\,=-\,\left(\frac{{k}_{\mathrm{hnf}}}{{k}_{f}}+\frac{4}{3}R\right)\left(\frac{1}{1-{s}_{t}}\right)\theta ^{\prime} \left(0\right),\end{array}\end{eqnarray}$
where $\mathrm{Re}=\tfrac{xu}{{v}_{f}}$ is the local Reynolds number.
3. Numerical method
The governing set of ordinary differential equations (10) and (11 ) along with boundary conditions (12) is numerically integrated by employing the built-in MATLAB solver bvp4c. The flow and heat transport characteristics are observed by the virtue of non-dimensional velocity, temperature, skin-friction, and heat transfer rate computed by the above-mentioned numerical method. The main purpose of this analysis is to predict multiple solutions for flow fields in the case of a shrinking surface. These solutions can be achieved by two different sets of initial guesses which fulfil the far-field boundary conditions. In this method, higher-order non-linear ODEs are first converted into first-order differential equations by choosing suitable variables. For this, let us assume $f={y}_{1},\,$ $f^{\prime} ={y}_{2},\,$ $f^{\prime\prime} ={y}_{3},\,$ $\theta ={y}_{4}\,$ and $\theta ^{\prime} ={y}_{5}\,.$. The new system of first-order differential equations is given by
$\begin{eqnarray}\begin{array}{l}{y^{\prime} }_{1}={y}_{2},\,{y^{\prime} }_{2}={y}_{3},\,\\ {y^{\prime} }_{3}=f^{\prime\prime\prime} ={(1-{\phi }_{1})}^{2.5}{(1-{\phi }_{2})}^{2.5}\\ \,\times \left\{(1-{\phi }_{2})\left[1-{\phi }_{1}+{\phi }_{1}\left(\displaystyle \frac{{\rho }_{s1}}{{\rho }_{f}}\right)\right]+{\phi }_{2}\left(\displaystyle \frac{{\rho }_{s2}}{{\rho }_{f}}\right)\right\}\,\\ \,\times \left[{y}_{2}^{2}-{y}_{1}{y}_{3}-1+\beta \left({y}_{2}+\displaystyle \frac{1}{2}\eta {y}_{3}-1\right)\right]\\ \,-M{(1-{\phi }_{1})}^{2.5}{(1-{\phi }_{2})}^{2.5}\left(1-{y}_{2}\right),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{y}^{{\prime} }_{4}={y}_{5},\,\\{y}^{{\prime} }_{5}=\theta ^{\prime\prime} =-\frac{\Pr }{\left(\frac{{k}_{\mathrm{hnf}}}{{k}_{f}}+\frac{4}{3}R\right)}\\ \,\times \left\{\left(1-{\phi }_{2}\right)\,\left[1-{\phi }_{1}+{\phi }_{1}\frac{{\left(\rho {C}_{p}\right)}_{s1}}{{\left(\rho {C}_{p}\right)}_{f}}\right]+{\phi }_{2}\frac{{\left(\rho {C}_{p}\right)}_{s2}}{{\left(\rho {C}_{p}\right)}_{f}}\right\}\,\\ \,\times \left[{y}_{1}{y}_{5}-{y}_{2}\left({y}_{4}+{s}_{t}\right)-\beta \left({y}_{2}+\frac{\eta }{2}{y}_{5}+{s}_{t}\right)\right]+\delta {y}_{4},\end{array}\end{eqnarray}$
and the boundary conditions become $\begin{eqnarray}\begin{array}{l}{y}_{1}\left(0\right)=S,\,{y}_{2}\left(0\right)=\lambda {S}_{1}-{\alpha }_{1}{y}_{3}\left(0\right),\,{y}_{4}\left(0\right)\\ \,\,\,\,=\,1-{s}_{t}\,,{y}_{2}\left(\infty \right)\to 1\,{\rm{and}}\,{y}_{4}\left(\infty \right)\to 0.\end{array}\end{eqnarray}$
An important step in this methodology is to give a suitable finite value to $(\eta =\infty )\,,$ so that the field boundary conditions (17) are asymptotically satisfied. The tolerance criterion is set as ${10}^{-6}\,$ to gain the correct numerical outcomes.4. Computed results
4.1. Code validation
To begin, tables 2 and 3 highlight the validity of the current results. In these table, the computed results of the skin-friction coefficient (upper and lower branch solutions) are compared with the published works of Wang et al [39], Mahapatra et al [40], Ismail et al [36] and Bachok et al [41]. In this regard, remarkable consistency has been accomplished. This proves that the proposed model and present findings are accurate.
Table 2. Validation of the computed results $f^{\prime\prime} (0)$ with existing works for distinct values of $\lambda $ when $S=M=0={\phi }_{1}={\phi }_{2}={\alpha }_{1}=\beta .$ |
| Wang et al [39] | Mahapatra et al [40] | Ismail et al [36] | Present study | |||||
|---|---|---|---|---|---|---|---|---|
| $\lambda $ | 1st sol. | 2nd sol. | 1st sol. | 2nd sol. | 1st sol. | 2nd sol. | 1st sol. | 2nd sol. |
| −0.25 | 1.402 24 | — | 1.402 242 | — | 1.402 240 7 | — | 1.402 240 189 17 | — |
| −1.00 | 1.328 82 | 0.0 | 1.328 819 | 0.0 | 1.328 816 8 | 0.0 | 1.328 817 063 98 | 0.0 |
| −1.15 | 1.082 23 | 0.116 702 | 1.082 232 | 0.116 702 | 1.082 231 1 | 0.116 672 4 | 1.082 231 659 85 | 0.116 702 098 941 |
| −1.20 | — | — | 0.932 470 | 0.233 648 | 0.932 473 3 | 0.233 628 4 | 0.932 473 336 374 | 0.233 649 713 559 |
Table 3. Comparison of results for upper and lower branch solution of $f^{\prime\prime} (0)$ for copper–water nanofluid $\left({\phi }_{1}=0\right)$ with variation in ${\phi }_{2}$ and $\lambda $ when $S=0=M=\,\,{\alpha }_{1}.$ |
| Bachok et al [41] | Present study | Bachok et al [41] | Present study | |||||
|---|---|---|---|---|---|---|---|---|
| $\lambda $ | ${\phi }_{2}=0$ | ${\phi }_{2}=0$ | ${\phi }_{2}=0.1$ | ${\phi }_{2}=0.1$ | ||||
| 1st sol. | 2nd sol. | 1st sol. | 2nd sol. | 1st sol. | 2nd sol. | 1st sol. | 2nd sol. | |
| 2 | −1.887 307 | — | −1.887 306 042 99 | — | −2.217 106 | — | −2.217 104 650 18 | — |
| 1 | 0 | — | 0 | — | 0 | — | 0 | — |
| 0.5 | 0.713 295 | — | 0.713 294 691 244 | — | 0.837 940 | — | 0.837 939 760 888 | — |
| 0 | 1.232 588 | — | 1.232 587 137 01 | — | 1.447 977 | — | 1.447 976 378 46 | — |
| −0.5 | 1.495 670 | — | 1.495 669 414 89 | — | 1.757 032 | — | 1.757 031 241 15 | — |
| −1 | 1.328 817 | 0 | 1.328 817 063 98 | — | 1.561 022 | 0 | 1.561 022 618 78 | — |
| −1.15 | 1.082 231 | 0.116 702 | 1.082 231 659 85 | 0.116 702 098 941 | 1.271 347 | 0.137 095 | 1.271 347 637 5 | 0.137 095 345 211 |
| −1.2 | 0.932 473 | 0.233 650 | 0.932 473 336 374 | 0.233 649 713 559 | 1.095 419 | 0.274 479 | 1.095 422 329 9 | 0.274 478 788 987 |
| −1.2465 | 0.584 281 | 0.554 297 | 0.584 282 077 652 | 0.554 296 057 004 | 0.686 379 | 0.651 161 | 0.686 382 925 642 | 0.651 157 307 193 |
4.2. Discussion
Various numerical and graphical outcomes of the problem have been prepared and presented in detail. The simulated outcomes are displayed through velocity and temperature distributions together with the skin-friction coefficient and heat transfer rate for different flow parameters to better understand the flow and heat transport features, for instance, for magnetic parameter $M,\,$ stretching/shrinking parameter $\lambda \,,$ suction parameter $S\,,$ velocity slip parameter ${\alpha }_{1},\,$ nanoparticle volume fractions ${\phi }_{1}\,$ and ${\phi }_{2}\,,$ unsteadiness parameter $\beta ,$ radiation parameter $R,\,$ thermal stratification parameter ${s}_{t},$ heat generation/absorption parameter $\delta $ and Prandtl number $\Pr .$
Figure 2 is plotted to witness the impact of magnetic parameter $M\,$ on wall shear stress $f^{\prime\prime} (0)\,$ against shrinking parameter $\lambda .$ It is shown that a dual solution exists for shrinking parameter $\left(\lambda \lt -1\right)\,$ and suction parameter $\left(S\gt 2\right)\,$ when keeping all other parameters fixed. Further, we observed that these dual solutions occur up to a critical value ${\lambda }_{c}\,$ of the shrinking parameter $\lambda .\,$ The critical value ${\lambda }_{c}\,$ plays an important role in predicting the nature of computed solutions. It is seen that a unique solution exists when $\lambda ={\lambda }_{c}\,$ and dual solutions are possible for $\lambda \lt {\lambda }_{c}\,$ while no solution exists when $\lambda \gt $ ${\lambda }_{c}\,.$ The respective critical values ${\lambda }_{c}\,$ for distinct magnetic parameters $M=(0,0.2,0.6)$ are ${\lambda }_{c}=-3.3556,\,$ $-3.3840\,$ and $-3.4019,\,$ as displayed in figure 2. A significant decline in the critical value is noted with increasing value of the magnetic parameter which enhances the existence domain of dual solutions. From a physical point of view, this happens because greater values of the magnetic parameter create resistance in the flow, which results in smothering of the velocity, as shown in the first solution. Moreover, the value of $f^{\prime\prime} (0)\,$ with respect to the first solution is always increasing with increasing magnetic parameter. On the other hand, an opposite behaviour is depicted by the second solutions. The impact of radiation parameter $R\,$ on heat transfer rate $-\theta ^{\prime} (0)\,$ against shrinking parameter $\lambda \,$ is displayed in figure 3. Again, we depict the existence of dual solutions in the case of shrinking flow. The magnitude of the critical value $\left|{\lambda }_{c}\right|\,$ increases with increasing radiation parameter $R.\,$ The wall shear stress $f^{\prime\prime} (0)\,$ is plotted in figure 4 for distinct values of suction parameter $S\,$ as a function of shrinking parameter $\lambda .\,$ It is worth mentioning that dual solutions exists for the shrinking parameter $(\lambda \lt 0)\,$ within a specific range of the suction parameter, i.e. $S=2.1,\,2.3,\,2.5\,.$ One can clearly see that the local wall shear stress $f^{\prime\prime} (0)\,$ increases with increasing suction parameter for the first solution, while it decreases in the case of the second solution. Figure 5 portrays the variation of local heat transfer rate $-\theta ^{\prime} (0)\,$ with varying values of thermal stratification parameter ${s}_{t}\,.$ We observe that as the stratification parameter increases the rate of heat transfer decreases for both solutions.
Figure 2. Variation of $f^{\prime\prime} (0)\,$ with magnetic parameter $M\,$ against $\lambda .$ |
Figure 3. Variation of $-\theta ^{\prime} \left(0\right)\,$ with radiation parameter $R\,$ against $\lambda .$ |
Figure 4. Variation of $f^{\prime\prime} (0)\,$ with suction parameter $S\,$ against $\lambda .$ |
Figure 5. Variation of $-\theta ^{\prime} \left(0\right)\,$ with stratification parameter ${s}_{t}\,$ against $\lambda .$ |
The impacts of the volume fraction of copper nanoparticles ${\phi }_{1}\,$ on the velocity and temperature distributions for fixed values of other parameter are displayed in figures 6 and 7. The outcomes uncovered that the variation in nanoparticle volume fraction has less effect on velocity fields. However, figure 6 suggests that the fluid velocity shows a decreasing behaviour with higher nanoparticle volume fraction for both solutions. On the other hand, it is observed through figure 7 that a substantial rise in temperature distribution is noted for greater volume fraction. The velocity distributions of $f^{\prime} (\eta )\,$ for various values of suction parameter $S\,$ are plotted in figure 8. The plotted graphs show the existence of dual solutions for shrinking flow inside the boundary layer region. As expected, an enhancement in velocity profiles is seen for greater $S\,$ in the upper solutions, while the inverse is noted in the lower solutions. Figure 9 portrays the relationship between the dimensionless temperature profiles $\theta (\eta )\,$ and the nanoparticle volume fraction ${\phi }_{2}\,$ within the boundary layer. It should be noted that both the temperature profiles and associated boundary layer thickness increase with growing values of ${\phi }_{2}.\,$ The velocity distributions $f^{\prime} (\eta )\,$ for varying $M\,$ are illustrated in figure 10, which shows that as the magnetic parameter $M\,$ increases the velocity profiles show a decreasing trend for the second solution and the opposite is noted for the first solutions. The dual temperature profiles $\theta (\eta )\,$ for precise entries of the radiation parameter $R\,$ are presented in figure 11. With an increment in radiation parameter, both solutions showed that the temperature distribution increases. On the other hand, the corresponding thermal boundary layer thickness is higher for the second solution. Figure 12 depicts the effect of velocity-slip parameter ${\alpha }_{1}\,$ on the dual velocity profiles $f^{\prime} (\eta )\,$ inside the boundary layer region. Once again, an opposite behaviour is displayed by the first and second branch solutions within the boundary layer region. The larger values of ${\alpha }_{1}\,$ cause a reduction in fluid velocity for the first branch and enhance the fluid velocity in the second branch. Finally, figure 13 is outlined for the behaviour temperature distribution $\theta (\eta )\,$ under the influence of stratification parameter ${s}_{t}\,.$ The hybrid nanofluid temperature is seen to rise for increasing values of the thermal slip parameter.
Figure 6. Velocity fields $f^{\prime} \left(\eta \right)$ for distinct nanoparticle volume fraction ${\phi }_{1}.$ |
Figure 7. Temperature fields $\theta \left(\eta \right)\,$ for distinct nanoparticle volume fraction ${\phi }_{1}.$ |
Figure 8. Velocity fields $f^{\prime} \left(\eta \right)$ for distinct suction parameter $S\,.$ |
Figure 9. Temperature fields $\theta \left(\eta \right)\,$ for distinct nanoparticle volume fraction ${\phi }_{2}\,.$ |
Figure 10. Velocity fields $f^{\prime} \left(\eta \right)$ for distinct magnetic parameter $M\,.$ |
Figure 11. Temperature fields $\theta \left(\eta \right)\,$ for distinct radiation parameter $R\,.$ |
Figure 12. Velocity fields $f^{\prime} \left(\eta \right)$ for distinct velocity-slip parameter ${\alpha }_{1}.$ |
Figure 13. Temperature fields $\theta \left(\eta \right)\,$ for distinct stratification parameter ${s}_{t}.$ |
5. Main findings
Numerical simulations for an unsteady thermally stratified flow of MHD hybrid nanofluid across a stretching/shrinking surface with thermal radiation and slip mechanism were carried out in this study. The main feature of the current study was the depiction of multiple branches, namely the upper and lower branch, for flow and temperature fields in the case of a shrinking surface. The following major outcomes can be summarized.
| 1. As the suction parameter was increased, the existence domain of the dual solution was increased with higher critical values of the shrinking parameter. | |
| 2. A higher skin-friction coefficient was noted for larger values of the magnetic parameter in the upper branch. | |
| 3. At higher values of the stratification parameter, a substantial rise in fluid temperature was observed for both solutions. | |
| 4. Higher values of the radiation and thermal stratification parameters decreased the Nusselt number for both the upper and lower branch solutions. | |
| 5. A decreasing tendency was observed for velocity curves with increased values of the velocity-slip parameter in the case of the second solution. | |
| 6. The hybrid nanofluid temperature was significantly increased by a greater thermal radiation parameter in both solutions. |