Radiative MHD Sutterby Nanofluid Flow Past a Moving Sheet: Scaling Group Analysis

Abstract

The present theoretical work endeavors to solve the Sutterby nanofluid flow and heat transfer problem over a permeable moving sheet, together with the presence of thermal radiation and magnetohydrodynamics (MHD). The fluid flow and heat transfer features near the stagnation region are considered. A new form of similarity transformations is introduced through scaling group analysis to simplify the governing boundary layer equations, which then eases the computational process in the MATLAB bvp4c function. The variation in the values of the governing parameters yields two different numerical solutions. One of the solutions is stable and physically reliable, while the other solution is unstable and is associated with flow separation. An increased effect of the thermal radiation improves the rate of convective heat transfer past the permeable shrinking sheet.

1. Introduction

The viscoelastic fluid is a type of non-Newtonian fluid that manifests the viscous and elasticity features under deformation. Sutterby fluid is an example of the viscoelastic fluid, and it well portrays the dilute polymer solutions [1,2]. Specifically, the Sutterby model fluid resembles the shear thinning and shear thickening aspects in high polymer aqueous solutions such as carboxymethyl cellulose (CMC), hydroxyethyl cellulose (HEC) and methyl cellulose (MC) [3]. The dilute polymer solutions have a wide range of functions in industrial practice, for instance, spray applications of agricultural chemicals [4], drag reducers in pipe flows [5], and production of domestic cleaning products [6]. The work of Fujii et al. [7] is one of the earliest studies to address the natural convection boundary layer flow in a Sutterby fluid past a vertical motionless isothermal plane and achieved an excellent comparison with the experimental results. Fujii et al. [8] revisited their work in [7] to investigate the impact of uniform heat flux under the same settings. However, the Sutterby model fluid received less attention from the boundary layer researchers at that time. Later, a new type of heat conductive fluid was introduced by Choi [9] named nanofluid. Nanofluid was also claimed to be a brilliant fluid due to its excellent heat transfer performance in engineering applications such as cooling of electronic appliances, and systems of solar water heating [10]. Nanofluid has now attracted significant interest from researchers, and boundary layer models have been studied under various settings [11,12,13]. After a long discontinuity in the theoretical works of the Sutterby boundary layer fluid flow, some numerical investigations of the Sutterby fluid under the Cattaneo–Christov heat flux [14,15], Soret and Dufour effect [16], peristaltic flow [17,18,19], squeezed flow [20], Joule heating effect [21], homogeneous–heterogeneous reactions [22], and hybrid nanoparticles [23] have been reported recently.

Magnetohydrodynamics (MHD) is another technological conception that is widespread in engineering practice. Electromagnetic casting [24], plasma confinement, and MHD power generation [25] are examples of notable applications. Thermal radiation is a type of energy that works in conjunction with the MHD effect. Thermal radiation emits and absorbs energy in the form of waves or molecules through a non-scattering medium. The successful combination of thermal radiation and MHD in an electrically conducting fluid has significant applications in solar power technology and electrical power generation [26]. Acknowledging these applications, researchers began to examine thermal radiation and MHD effects in the boundary layer flow past a stretching/shrinking surface, and many theoretical works have been reported. Recently, Sabir et al. [27] explored the stagnation-point flow of a Sutterby fluid with the effects of an inclined magnetic field and thermal radiation past a stretching surface, and observed the declining trend of the convective heat transfer with the stronger influence of thermal radiation. Bilal et al. [28] examined the ohmically dissipated Darcy–Forchheimer slip flow of an MHD Sutterby fluid past a radiating stretching sheet and found a decrement in the convective heat transfer with increasing slip effects.

By comparison, the boundary layer equations, which were proposed by Prandtl [29], disclosed many invariant closed-form solutions. Prandtl’s boundary layer equations can be reduced to a less complicated form that is in a system of ordinary differential equations. These boundary layer equations also allow many different types of symmetry groups, of which the Lie group analysis is prominent. Lie group analysis helps to identify the transformation point that represents the given boundary layer equations [30]. In Lie group analysis, the group-invariant solutions are the similarity solutions, and these similarity solutions are used to reduce the independent variables in a fluid flow problem [31]. A special form of the Lie group analysis exists, namely, the scaling group of transformation, and this has been employed by researchers in valuable contributions, for instance, see [32,33].

Regarding studies of stagnation-point flow in a Sutterby fluid, Azhar et al. [34] investigated the effect of entropy generation on the stagnation-point flow of a Sutterby nanofluid past a stretching sheet. Azhar et al. [35] reconsidered the work of [34] by incorporating the Cattaneo–Christov heat flux model and omitting the nanoparticles. Both of the studies of [34,35] solved the flow problem numerically and presented unique solutions. A number of considerable research gaps were found in the theoretical works available in the stagnation-point flow and heat transfer in a Sutterby nanofluid, for instance: inspecting fluid flow behavior and heat transfer characteristics past a shrinking sheet together with the suction effect; conducting scaling group analysis; obtaining dual solutions; and performing stability analysis. Thus, the present work is devoted to numerically solving the problem of boundary layer Sutterby nanofluid flow and heat transfer near the stagnation region over a permeable moving (stretching/shrinking) sheet. The fluid flow and heat transfer characteristics under the magnetic and thermal radiation effects are observed. Scaling group analysis is employed to obtain the apt similarity transformations so that the complex governing boundary layer equations can be brought to a soluble form. The simplified form of the mathematical model is then solved numerically in the boundary value problem solver or bvp4c function in MATLAB. Two different numerical solutions are identified with the governing parameters’ variation. Further, stability analysis is undertaken in the present work to justify the presence of dual solutions. These contributions are essentially original, and all numerical results are presented and discussed in detail.

2. Problem Formulation

Contemplate an incompressible two-dimensional stagnation-point flow of a Sutterby nanofluid across a stretching/shrinking sheet as shown in Figure 1, where x and y are the Cartesian coordinates with the x-axis positioned in the horizontal direction, and the y-coordinate is normal to the x-coordinate. The free stream velocity is denoted by ${\overline{u}}_e,$ and ${\overline{u}}_w$ signifies the velocity of the moving sheet, where ${\overline{u}}_w>0$ infers the state of the stretching sheet, ${\overline{u}}_w<0$ the embodies shrinking sheet, and ${\overline{u}}_w=0$ typifies the stationary sheet. The moving (stretching or shrinking) sheet is penetrable and there is a uniform surface mass flux, of velocity ${\overline{v}}_w,$ with ${\overline{v}}_w>0$ to imply the injection situation and ${\overline{v}}_w<0$ for the suction state. The free stream temperature and the wall temperature are denoted by $T_{\infty }$ and $T_w,$ respectively.

Sutterby [1,2] introduced the constitutive law for the Sutterby fluid by expressing the Cauchy stress tensor $\left(\mathbf{T}\right)$ as:

$$\mathbf{T}=-p\mathbf{I}+\mathbf{S},$$

(1)

where $p$ is the pressure, $\mathbf{I}$ is the identity vector, and $\mathbf{S}$ is the extra stress tensor which can be defined as follows [21]:

$$\mathbf{S}={\mu }_0{\left[\frac{{\sinh }^{-1}\left(E\dot{\gamma }\right)}{\left(E\dot{\gamma }\right)}\right]}^m{\mathbf{A}}_1.$$

(2)

Here, ${\mu }_0$ is the viscosity at low shear rates, $E$ is the material time constant, $\dot{\gamma }=\sqrt{\mathrm{tr}{\left({\mathbf{A}}_1\right)}^2/2}$ is the second invariant strain tensor, ${\mathbf{A}}_1$ is the first order Rivlin–Erickson tensor or deformation rate tensor which is defined as ${\mathbf{A}}_1=\left(\nabla \mathbf{V}\right)+{\left(\nabla \mathbf{V}\right)}^{\mathrm{T}},$ and $m$ is the power-law index. The Sutterby model in Equation (2) is a versatile model when the value of $m$ changes. For instance, when $m=0,$ the Sutterby model imitates the Newtonian fluid behavior, when $m=1,$ the model is reduced to the Eyring model, and this model also predicts specifically the pseudo-plastic (shear thinning) and dilatant (shear thickening) fluid properties when $m>0$ and $m<0,$ respectively. By taking the velocity field as $\mathbf{V}=\left[u\left(x,y\right),v\left(x,y\right)\right],$ and under the assumptions mentioned earlier, the governing boundary layer equations in the dimensional form can be formed as follows [36]:

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

(3)

$$\overline{u}\frac{\partial \overline{u}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{u}}{\partial \overline{y}}=\frac{{\mu }_{nf}}{{\rho }_{nf}}\frac{{\partial }^2\overline{u}}{\partial {\overline{y}}^2}+\frac{{\mu }_{nf}}{{\rho }_{nf}}\frac{m{E}^2}{2}{\left(\frac{\partial \overline{u}}{\partial \overline{y}}\right)}^2\frac{{\partial }^2\overline{u}}{\partial {\overline{y}}^2}+{\overline{u}}_e\frac{d{\overline{u}}_e}{d\overline{x}}-\frac{\sigma B_0^2\left({\overline{u}}_e-\overline{u}\right)}{{\rho }_{nf}},$$

(4)

$$\overline{u}\frac{\partial T}{\partial \overline{x}}+\overline{v}\frac{\partial T}{\partial \overline{y}}=\frac{k_{nf}}{{\left(\rho C_p\right)}_{nf}}\frac{{\partial }^{2}T}{\partial {\overline{y}}^2}+\frac{16{\sigma }_1{T}_{\infty }^3}{3{\left(\rho C_p\right)}_{nf}{k}_1}\frac{{\partial }^{2}T}{\partial {\overline{y}}^2},$$

(5)

along with the respective boundary conditions:

$$\begin{array}{l}\overline{u}={\overline{u}}_w,\overline{v}={\overline{v}}_w,T=T_w\left(x\right)\mathrm{at}\overline{y}=0.\\ \overline{u}={\overline{u}}_e,\frac{\partial \overline{u}}{\partial \overline{y}}\to 0,T\to T_{\infty }\mathrm{as}\overline{y}\to \infty ,\end{array}$$

(6)

where $\overline{u}$ and $\overline{v}$ denote velocity components in the $\overline{x}$ and $\overline{y}$ directions, respectively, ${\mu }_{nf}$ is the dynamic viscosity of the nanofluid, $\sigma$ is the electrical conductivity, $B_0$ is the magnetic field strength, ${\rho }_{nf}$ is the density of the nanofluid, ${\sigma }_1$ is the Stefan Boltzmann constant, $k_1$ is the Rosseland mean absorption coefficient, $k_{nf}$ is the thermal conductivity of the nanofluid, and ${\left(C_p\right)}_{nf}$ is the specific heat capacity of the nanofluid. The detailed definitions of the nanofluid parameters are given by the following expressions, which are valid when the nanoparticles are of spherical shape or similar to a spherical shape [37]:

$$\begin{array}{l}{\mu }_{nf}=\frac{{\mu }_{bf}}{{\left(1-\varphi \right)}^{2.5}},{\alpha }_{nf}=\frac{k_{nf}}{{\left(\rho C_p\right)}_{nf}},{\left(\rho C_p\right)}_{nf}=\left(1-\varphi \right){\left(\rho C_p\right)}_{bf}+\varphi {\left(\rho C_p\right)}_s,\\ {\rho }_{nf}=\left(1-\varphi \right){\rho }_{bf}+\varphi {\rho }_s,\frac{k_{nf}}{k_{bf}}=\frac{\left(k_s+2k_{bf}\right)-2\varphi \left(k_{bf}-k_s\right)}{\left(k_s+2k_{bf}\right)+\varphi \left(k_{bf}-k_s\right)},\end{array}$$

(7)

where $\varphi$ denotes the nanoparticle volume fraction, ${\mu }_{bf}$ denotes the dynamic viscosity of the base fluid, ${\alpha }_{nf}$ denotes the thermal diffusivity of the nanofluid, $k_{bf}$ is the thermal conductivity of the base fluid, $k_s$ is the thermal conductivity of the solid fractions, $C_p$ is the specific heat capacity, and ${\rho }_{bf}$ and ${\rho }_s$ are the densities of the base fluid and solid fractions, respectively. The Sutterby model reflects the dilute polymer solution where the polymer is diluted in the appropriate solvent. Hence, for the present study, n-Hexane is chosen as the base fluid (solvent).

Table 1. The thermo physical characteristics of the essential fluid and nanoparticles.

Physical Properties	Fluid Phase (n-Hexane /CH3(CH2)4CH3)	Solid Phase Magnetite (Fe3O4)
$C_p\left(\mathrm{kJ}/\mathrm{kg}\cdot \mathrm{K}\right)$	2.78	670
$\rho \left(\mathrm{kg}/{\mathrm{m}}^3\right)$	551	5180
$k\left(\mathrm{W}/\mathrm{mK}\right)$	82	9.7
Pr	4.36	$-$

displays the specific values for the respective thermophysical features of n-Hexane and magnetite nanofluid [38].

3. Non-Dimensionalization of the Governing Equations

Considering the following the non-dimensional variables:

$$\begin{array}{l}x=\frac{a\overline{x}}{u_0},y=\sqrt{\frac{a}{{\nu }_{bf}}}\overline{y},u=\frac{\overline{u}}{u_0},u_w=\frac{{\overline{u}}_w}{u_0},v=\frac{\overline{v}}{\sqrt{a{\nu }_{bf}}},\\ v_w=\frac{{\overline{v}}_w}{\sqrt{a{\nu }_{bf}}},u_e=\frac{{\overline{u}}_e}{u_0},\theta =\frac{T-T_{\infty }}{T_0},\end{array}$$

(8)

where $u_0$ is the characteristic velocity and introducing the stream function $\psi ,$ which can be defined by $u=\frac{\partial \psi }{\partial y}$ and $v=-\frac{\partial \psi }{\partial x}$, Equations (4) and (5) become:

$$\frac{\partial \psi }{\partial y}\frac{{\partial }^2\psi }{\partial x\partial y}-\frac{\partial \psi }{\partial x}\frac{{\partial }^2\psi }{\partial y^2}=\frac{A_1}{A_2}\frac{{\partial }^3\psi }{\partial y^3}+\frac{A_1}{A_2}\frac{mDe}{2}{\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}^2\frac{{\partial }^3\psi }{\partial y^3}+u_e\frac{d{u}_e}{dx}-\frac{\sigma B_o^2}{{\rho }_{bf}a}\frac{1}{A_2}\left(u_e-\frac{\partial \psi }{\partial y}\right),$$

(9)

$$\frac{\partial \psi }{\partial y}\frac{\partial \theta }{\partial x}-\frac{\partial \psi }{\partial x}\frac{\partial \theta }{\partial y}=\frac{A_4}{A_3}\frac{1}{Pr}\frac{{\partial }^2\theta }{\partial y^2}+\frac{1}{A_3}\frac{4}{3}\frac{Rd}{Pr}\frac{{\partial }^2\theta }{\partial y^2},$$

(10)

with the corresponding boundary conditions:

$$\begin{array}{l}\frac{\partial \psi }{\partial y}=u_w,\frac{\partial \psi }{\partial x}=-v_w,\theta T_0=T_w\left(x\right)\mathrm{at}y=0,\\ \frac{\partial \psi }{\partial y}\to u_e,\frac{{\partial }^2\psi }{\partial y^2}\to 0,\theta \to 0\mathrm{as}y\to \infty ,\end{array}$$

(11)

while satisfying the continuity equation of Equation (3). In Equations (9) and (11), $M=\frac{\sigma B_0^2}{{\rho }_{bf}a}$ is the magnetic parameter, $Rd=\frac{4{\sigma }_1{T}_{\infty }^3}{k_1{k}_{bf}}$ is the radiation parameter, $\mathrm{Pr}=\frac{{\mu }_{bf}{\left(c_p\right)}_{bf}}{k_{bf}}$ is the Prandtl number, $De=\frac{u_0^{2}a{E}^2}{{\nu }_{bf}}$ is the Deborah number, $\varphi$ is the nanoparticle volume fraction, and terms $A_1,A_2,A_3,$ and $A_4$ are expressed as:

$$\begin{array}{l}{A}_1=\frac{1}{{\left(1-\varphi \right)}^{2.5}},A_2=1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_{bf}},A_3=1-\varphi +\varphi \frac{{\left(\rho c_p\right)}_s}{{\left(\rho c_p\right)}_{bf}},\\ A_4=\frac{k_{nf}}{k_f}=\frac{k_s+2k_{bf}-2\varphi \left(k_{bf}-k_s\right)}{k_{bf}+2k_{bf}+\varphi \left(k_{bf}-k_s\right)}.\end{array}$$

(12)

The functions $u_w,v_w$ and $T_w\left(x\right)$ are assumed to be in the following form to ensure that similarity solution exists:

$$u_w=\frac{u_1}{u_0}{x}^{\frac{2}{5}},v_w=\frac{v_1}{\sqrt{a{\nu }_{bf}}}{x}^{-\frac{2}{5}},T_w\left(x\right)=T_0{x}^{\frac{2}{5}},$$

(13)

where $u_1$ is the reference velocity, $v_1$ is the normal reference velocity, and $T_0$ is the reference temperature.

4. Scaling Group Analysis

The governing boundary layer flow and heat transfer problem in the form of partial differential equations (PDEs) is complex and hard to solve by means of mathematical software. Therefore, it needs to be reduced to a simpler form so that it can be solved. Suitable similarity variables can facilitate the transformation and, at this point, scaling group analysis is required to form the specified similarity transformations for the present problem. The newly formed similarity variable will then transform the PDEs to a system of ordinary differential equations (ODEs), and the model can be solved by the desired mathematical software. Therefore, the following scaling group of transformations G is introduced:

$$\begin{array}{l}G:x^{*}=x{G}^{{\omega }_1},y^{*}=y{G}^{{\omega }_2},{\psi }^{*}=\psi G^{{\omega }_3},{\sigma }^{*}=\sigma G^{{\omega }_4},\\ {\theta }^{*}=\theta G^{{\omega }_5},u_e^{*}=u_e{G}^{{\omega }_6},u_1^{*}=u_1{G}^{{\omega }_7},m^{*}=m{G}^{{\omega }_8},\end{array}$$

(14)

where ${\omega }_i$ are constants to be determined in which $i=1,\dots 8.$ The transformation $G$ is the transformation point which transforms the $\left(x,y,\psi ,\sigma ,\theta ,u_e,u_1,m,\right)$ coordinates to the new coordinates $\left(x^{*},y^{*},{\psi }^{*},{\sigma }^{*},{\theta }^{*},u_e^{*},u_1^{*},m^{*}\right).$

Next, the substitution of (14) into Equations (9)–(11) yields the following expressions:

$$\begin{array}{l}\frac{A_1}{A_2}{G}^{\left[2{\omega }_3-2{\omega }_2-{\omega }_1\right]}\left(\frac{\partial {\psi }^{*}}{\partial y^{*}}\frac{{\partial }^2{\psi }^{*}}{\partial x^{*}\partial y^{*}}-\frac{\partial {\psi }^{*}}{\partial x^{*}}\frac{{\partial }^2{\psi }^{*}}{\partial y^{*2}}\right)=G^{\left[{\omega }_3-3{\omega }_2\right]}\left(\frac{{\partial }^3{\psi }^{*}}{\partial y^{*3}}\right)\\ +\frac{A_1}{A_2}\frac{De}{2}{G}^{\left[{\omega }_8+3{\omega }_3-7{\omega }_2\right]}\left[m^{*}{\left(\frac{{\partial }^2{\psi }^{*}}{\partial y^{*2}}\right)}^2\frac{{\partial }^3{\psi }^{*}}{\partial y^{*3}}\right]+G^{\left[2{\omega }_6-{\omega }_1\right]}{u}_e^{*}\frac{d{u}_e^{*}}{d{x}^{*}}\\ -\frac{B_0^2}{{\rho }_{bf}a}\frac{1}{A_2}\left[G^{\left[{\omega }_3+{\omega }_4-{\omega }_2\right]}\left({\sigma }^{*}\frac{\partial {\psi }^{*}}{\partial y^{*}}\right)-G^{\left[{\omega }_4+{\omega }_6\right]}\left({\sigma }^{*}{u}_e^{*}\right)\right],\end{array}$$

(15)

$$G^{\left[{\omega }_3+{\omega }_5-{\omega }_1-{\omega }_2\right]}\left(\frac{\partial {\psi }^{*}}{\partial y^{*}}\frac{\partial {\theta }^{*}}{\partial x^{*}}-\frac{\partial {\psi }^{*}}{\partial x^{*}}\frac{\partial {\theta }^{*}}{\partial y^{*}}\right)=\frac{1}{Pr}{G}^{\left[{\omega }_5-2{\omega }_2\right]}\left(\frac{{\partial }^2{\theta }^{*}}{\partial y^{*2}}\right)+\frac{4}{3}\frac{Rd}{Pr}{G}^{\left[{\omega }_5-2{\omega }_2\right]}\left(\frac{{\partial }^2{\theta }^{*}}{\partial y^{*2}}\right),$$

(16)

along with the boundary conditions:

$$\begin{array}{l}{G}^{\left[{\omega }_3-{\omega }_2\right]}\left(\frac{\partial {\psi }^{*}}{\partial y^{*}}\right)=G^{\left[{\omega }_7+\frac{2}{5}{\omega }_1\right]}\left(\frac{u_1^{*}}{u_0}{x}^{*\frac{2}{5}}\right),\\ G^{\left[{\omega }_3-{\omega }_1\right]}\left(\frac{\partial {\psi }^{*}}{\partial x^{*}}\right)=-\frac{v_1}{\sqrt{a{\nu }_{bf}}}{G}^{\left[-\frac{2}{5}{\omega }_1\right]}{x}^{*-\frac{2}{5}},G^{\left[{\omega }_5\right]}{\theta }^{*}=G^{\left[\frac{2}{5}{\omega }_1\right]}\left(x^{*\frac{2}{5}}\right)\mathrm{at}y=0\\ G^{\left[{\omega }_3-{\omega }_2\right]}\left(\frac{\partial {\psi }^{*}}{\partial y^{*}}\right)\to G^{\left[{\omega }_6\right]}\left(u_e^{*}\right)\mathrm{as}y\to \infty .\end{array}$$

(17)

To retain the invariance of the system under $G,$ the parameters defined in Equation (14), the following relations must hold:

$$\begin{array}{l}2{\omega }_3-2{\omega }_2-{\omega }_1={\omega }_3-3{\omega }_2=3{\omega }_3-7{\omega }_2+{\omega }_8=2{\omega }_6-{\omega }_1={\omega }_3-{\omega }_2+{\omega }_4={\omega }_4+{\omega }_6\\ ={\omega }_3+{\omega }_5-{\omega }_1-{\omega }_2={\omega }_5-2{\omega }_2.\end{array}$$

(18)

From the boundary conditions of Equations (17), we also obtain the following relations among the parameters:

$${\omega }_3-{\omega }_2={\omega }_7+\frac{2}{5}{\omega }_1,{\omega }_3-{\omega }_1=-\frac{2}{5}{\omega }_1,{\omega }_5=\frac{2}{5}{\omega }_1,{\omega }_3-{\omega }_2={\omega }_6.$$

(19)

The absolute invariant can be determined by eliminating the parameter G of the group and hence Equations (18) and (19) provide the following expressions:

$$\begin{array}{l}{\omega }_2=\frac{2}{5}{\omega }_1,{\omega }_3=\frac{3}{5}{\omega }_1,{\omega }_4=-\frac{4}{5}{\omega }_1,{\omega }_5=\frac{2}{5}{\omega }_1,\\ {\omega }_6=\frac{1}{5}{\omega }_1,{\omega }_7=-\frac{1}{5}{\omega }_1,{\omega }_8=\frac{2}{5}{\omega }_1.\end{array}$$

(20)

From Equations (13), (14), and (20), we achieve the absolute invariants under the group G similarity transformations as follows:

$$\begin{array}{l}\eta =\frac{y}{x^{\frac{2}{5}}},\psi =x^{\frac{3}{5}}f\left(\eta \right),\sigma ={\sigma }_0{x}^{-\frac{4}{5}},\theta ={\theta }_0\left(\eta \right)x^{\frac{2}{5}},\\ u_e={\left(u_e\right)}_0{x}^{\frac{1}{5}},u_1={\left(u_1\right)}_0{x}^{-\frac{1}{5}},m=m_0{x}^{\frac{2}{5}}.\end{array}$$

(21)

The similarity transformations in Equation (21) are new and, by employing them in the governing boundary layer equations of Equations (9) and (11), the reduced version of the model in the form of ordinary differential equations can be attained as follows while satisfying Equation (9):

$$\frac{A_1}{A_2}{f}^{‴}\left[1+\frac{m_{0}De}{2}{\left(f^{″}\right)}^2\right]-\frac{1}{5}{\left(f^{\prime }\right)}^2+\frac{3}{5}f{f}^{″}+\frac{1}{5}-\frac{M}{A_2}\left(1-f^{\prime }\right)=0,$$

(22)

$$\left(A_4+\frac{4}{3}Rd\right){\theta }^{″}-\frac{2}{5}{A}_3\mathrm{Pr}{f}^{\prime }\theta +\frac{3}{5}{A}_3\mathrm{Pr}f{\theta }^{\prime }=0,$$

(23)

with the associated boundary conditions:

$$f(0)=\frac{5}{3}{f}_w,f^{\prime }(0)=\epsilon ,\theta \left(0\right)=1,f^{\prime }(\infty )=1,\theta (\infty )=0.$$

(24)

Here $\epsilon ={\left(u_1\right)}_0/u_0$ is the stretching/shrinking parameter, where $\epsilon >0$ indicates the stretching sheet, $\epsilon =0$ specifies the stationary sheet, and $\epsilon <0$ represents the state of shrinking sheet. Furthermore, $f_w=-v_1/\sqrt{{\nu }_{bf}a}$ is the constant mass transfer parameter, and $f_w>0$ typifies the suction effect at the surface of the moving sheet and $f_w<0$ epitomizes the injection state. For simplicity, we choose ${\left(u_e\right)}_0=1.$ The power-law index is denoted by $m_0$; when $m_0=0,$ the reduced model imitates the Newtonian fluid behavior. Moreover, when $m_0=1,$ the model is reduced to the Eyring model, whereas the present model in Equations (22)–(24) also predicts specifically the pseudo-plastic (shear thinning) and dilatant (shear thickening) fluid properties when $m_0<0$ and $m_0>0,$ respectively.

The physical quantities of interest in the present work are the local skin friction coefficient $\left(C_{f\overline{x}}\right)$ and the local Nusselt number $\left(N{u}_{\overline{x}}\right)$ which are defined as follows:

$$C_{f\overline{x}}=\frac{{\tau }_w}{{\rho }_{bf}{\overline{u}}_e^2},N{u}_{\overline{x}}=\frac{\overline{x}{q}_w}{k_{bf}\left(T_w-T_{\infty }\right)},$$

(25)

where ${\tau }_w$ is the wall shear stress and $q_w$ is the heat flux at the surface of the sheet, and can be further defined as [34]:

$${\tau }_w={\mu }_{nf}\left\{1+\frac{m{E}^2}{6}\left[2{\left(\frac{\partial \overline{u}}{\partial \overline{x}}\right)}_{\overline{y}=0}^2+{\left(\frac{\partial \overline{u}}{\partial \overline{y}}\right)}_{\overline{y}=0}\right]\right\}{\left(\frac{\partial \overline{u}}{\partial \overline{y}}\right)}_{\overline{y}=0},q_w=-k_{nf}{\left(\frac{\partial T}{\partial \overline{y}}\right)}_{\overline{y}=0}.$$

(26)

The reduced skin friction coefficient ${\mathrm{Re}}_x^{1/2}{C}_{f\overline{x}}{x}^{1/5}$ and the local Nusselt number ${\mathrm{Re}}_x^{-1/2}N{u}_{\overline{x}}{x}^{-1/5}$ can be obtained using the similarity transformations of Equation (21) and the expressions in Equations (25) and (26) as follows:

$$\begin{array}{l}{C}_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}=\sqrt{A_1{A}_2}{f}^{″}(0)+\frac{m_0}{6\sqrt{A_1{A}_2}}De{\left[f^{″}(0)\right]}^3,\\ {\mathrm{Re}}_{\overline{x}}^{-1/2}N{u}_{\overline{x}}{x}^{-2/5}=-\frac{1}{\sqrt{A_1{A}_2}}\left(1+\frac{1}{A_4}\frac{4}{3}Rd\right)\frac{{\theta }^{\prime }\left(0\right)}{\theta \left(0\right)},\end{array}$$

(27)

where ${\mathrm{Re}}_{\overline{x}}={\overline{u}}_e\overline{x}/{\nu }_{bf}$ denotes the local Reynolds number.

5. Stability Analysis

Merkin [39,40] established an improved version of the stability analysis, which is prominent among researchers for the examination of the stability of numerical solutions. Because we observed dual solutions in the present work, we assess the solution’s stability to determine the flow behavior. To initiate the linear stability analysis, the model equations in Equations (3)–(5) need to be considered in the unsteady form as follows:

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

(28)

$$\frac{\partial \overline{u}}{\partial \overline{t}}+\overline{u}\frac{\partial \overline{u}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{u}}{\partial \overline{y}}=\frac{{\mu }_{nf}}{{\rho }_{nf}}\frac{{\partial }^2\overline{u}}{\partial {\overline{y}}^2}+\frac{{\mu }_{nf}}{{\rho }_{nf}}\frac{m{E}^2}{2}{\left(\frac{\partial \overline{u}}{\partial \overline{y}}\right)}^2\frac{{\partial }^2\overline{u}}{\partial {\overline{y}}^2}+{\overline{u}}_e\frac{d{\overline{u}}_e}{d\overline{x}}-\frac{\sigma B_0^2\left({\overline{u}}_e-\overline{u}\right)}{{\rho }_{nf}},$$

(29)

$$\frac{\partial T}{\partial \overline{t}}+\overline{u}\frac{\partial T}{\partial \overline{x}}+\overline{v}\frac{\partial T}{\partial \overline{y}}=\frac{k_{nf}}{{\left(\rho C_p\right)}_{nf}}\frac{{\partial }^{2}T}{\partial {\overline{y}}^2}+\frac{16{\sigma }_1{T}_{\infty }^3}{3{\left(\rho C_p\right)}_{nf}{k}_1}\frac{{\partial }^{2}T}{\partial {\overline{y}}^2},$$

(30)

with the boundary conditions of Equation (6). Then, we introduce the dimensional time variable, $\overline{t}$ with a new similarity variable $\left(\tau ={\tau }_0/x^{4/5}\right)$ through scaling group analysis. The new similarity transformation is given as:

$$\begin{array}{l}\eta =\frac{y}{x^{\frac{2}{5}}},\psi =x^{\frac{3}{5}}f\left(\eta ,\tau \right),\sigma ={\sigma }_0{x}^{-\frac{4}{5}},\theta ={\theta }_0\left(\eta ,\tau \right)x^{\frac{2}{5}},\\ u_e={\left(u_e\right)}_0{x}^{\frac{1}{5}},u_1={\left(u_1\right)}_0{x}^{-\frac{1}{5}},m=m_0{x}^{\frac{2}{5}},\tau =\frac{{\tau }_0}{x^{\frac{4}{5}}}.\end{array}$$

(31)

Employing (31) in the dimensionless form of Equations (28)–(30) and (6) gives the following system of equations:

$$\frac{A_1}{A_2}\frac{{\partial }^{3}f}{\partial {\eta }^3}\left[1+\frac{m_{0}De}{2}{\left(\frac{{\partial }^{2}f}{\partial {\eta }^2}\right)}^2\right]-\frac{1}{5}{\left(\frac{\partial f}{\partial \eta }\right)}^2+\frac{3}{5}f\left(\eta ,\tau \right)\frac{{\partial }^{2}f}{\partial {\eta }^2}+\frac{1}{5}-\frac{M}{A_2}\left(1-\frac{\partial f}{\partial \eta }\right)-\frac{{\partial }^{2}f}{\partial \eta \partial \tau }=0,$$

(32)

$$\left(A_4+\frac{4}{3}Rd\right)\frac{{\partial }^2\theta }{\partial {\eta }^2}-\frac{2}{5}{A}_3\mathrm{Pr}\frac{\partial f}{\partial \eta }\theta \left(\eta ,\tau \right)+\frac{3}{5}{A}_3\mathrm{Pr}f\left(\eta ,\tau \right)\frac{\partial \theta }{\partial \eta }-\frac{\partial \theta }{\partial \tau }=0,$$

(33)

with the boundary conditions:

$$f\left(0,\tau \right)=\frac{5}{3}{f}_w,\frac{\partial f}{\partial \eta }\left(0,\tau \right)=\epsilon ,\theta \left(0,\tau \right)=1,\frac{\partial f}{\partial \eta }\left(\infty ,\tau \right)=1,\theta \left(\infty ,\tau \right)=0.$$

(34)

It is assumed that the solutions of (32)–(34) are expressed by the formulas of Equation (35):

$$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),$$

(35)

where $f\left(\eta \right)=f_0\left(\eta \right)$ and $\theta \left(\eta \right)={\theta }_0\left(\eta \right)$ are the solutions found in the previous section, in which the disturbance is superimposed to determine their stability. Here, the unknown eigenvalue parameter is denoted by $\gamma$, and $F\left(\eta ,\tau \right)$ and $G\left(\eta ,\tau \right)$ are relatively small compared to the steady state solutions $(f_0\left(\eta \right)$ and ${\theta }_0\left(\eta \right)).$ The substitution of Equation (35) into Equations (32)–(34) gives the following system:

$$\frac{{\partial }^{3}F}{\partial {\eta }^3}\left(1+\frac{A_1}{A_2}\frac{m_{0}De}{2}{f_0^{″}}^2\right)+\frac{{\partial }^{2}F}{\partial {\eta }^2}\left(\frac{3}{5}{f}_0-\frac{2}{5}{f}_0^{\prime }\right)+\frac{3}{5}{f}_0^{″}F+\left(\frac{M}{A_2}+\gamma \right)\frac{\partial F}{\partial \eta }-\frac{{\partial }^{2}F}{\partial \eta \partial \tau }=0,$$

(36)

$$\left(A_4+\frac{4}{3}Rd\right)\frac{{\partial }^{2}G}{\partial {\eta }^2}+\frac{3}{5}{A}_3\mathrm{Pr}{f}_0\frac{\partial G}{\partial \eta }+\frac{3}{5}{A}_3\mathrm{Pr}F{\theta }_0-\frac{2}{5}{A}_3\mathrm{Pr}{f}_0^{\prime }G-\frac{2}{5}{A}_3\mathrm{Pr}{\theta }_0\frac{\partial F}{\partial \eta }+\gamma G-\frac{\partial G}{\partial \tau }=0,$$

(37)

subject to the boundary conditions:

$$F\left(0,\tau \right)=0,\frac{\partial F}{\partial \eta }\left(0,\tau \right)=0,G\left(0,\tau \right)=0,\frac{\partial F}{\partial \eta }\left(\infty ,\tau \right)=0,G\left(\infty ,\tau \right)=0.$$

(38)

Referring to Merkin [39,40], $\tau \to 0$ is fixed to examine the stability of the steady state boundary layer flow. Thus, $F=F_0\left(\eta \right)$ and $G=G_0\left(\eta \right)$ in Equations (37)–(39), yielding the following linearized eigenvalue problem:

$$F_0^{‴}\left(1+\frac{A_1}{A_2}\frac{m_{0}De}{2}{\left(f_0^{″}\right)}^2\right)+F_0^{″}\left(\frac{3}{5}{f}_0-\frac{2}{5}{f}_0^{\prime }\right)+\frac{3}{5}{f}_0^{″}{F}_0+\left(\frac{M}{A_2}+\gamma \right)F_0^{\prime }=0,$$

(39)

$$\left(A_4+\frac{4}{3}Rd\right)G_0^{″}+\frac{3}{5}{A}_3\mathrm{Pr}{f}_0{G}_0^{\prime }+\frac{3}{5}{A}_3\mathrm{Pr}{F}_0{\theta }_0-\frac{2}{5}{A}_3\mathrm{Pr}{f}_0^{\prime }{G}_0-\frac{2}{5}{A}_3\mathrm{Pr}{\theta }_0{F}_0^{\prime }+\gamma G_0=0,$$

(40)

with the boundary conditions:

$$\begin{array}{l}{F}_0\left(\eta \right)=0,F_0^{\prime }\left(\eta \right)=0,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}$$

(41)

It is necessary to replace one of the outer boundary conditions with a normalizing boundary condition to obtain the eigenvalues. Therefore, the boundary condition $F_0^{\prime }(\infty )=0$ is substituted with $F_0^{″}\left(0\right)=1.$ The system of equations in Equations (38)–(40) with the new boundary condition is solved by the MATLAB boundary value problem solver (bvp4c) to obtain the lowest eigenvalues as the governing parameter varies. These lowest eigenvalues are classified according to their sign. If the lowest eigenvalue falls in the positive range of values, then the respective numerical solution is accepted as a stable solution. Meanwhile, the negative lowest eigenvalue suggests the numerical solution is unstable. Further explanation about the stable and unstable solutions is provided in the next section.

6. Results and Discussion

The mathematical model in Equations (23)–(25) was solved numerically by means of the boundary value problem solver function bvp4c in the MATLAB software. The numerical results were derived while limiting the relative tolerance to $1\times {10}^{-10}.$ Some of the governing parameter values were fixed throughout the computation process to align with the motivation of this study. For example, the Prandtl number (Pr) value was fixed at 4.36 because it represents the base fluid, n-hexane. The power-law index $\left(m_0\right)$ was fixed at 1.5 to investigate the dilatant features of the Sutterby fluid. The obtained non-uniqueness solutions were classified based on how early the solution converged asymptotically. For example, the numerical solution that converged earlier asymptotically in the velocity/temperature profiles was labelled the first solution. The other solution that converged later, asymptotically, was labelled as the second solution. Before presenting the numerical results, we provide validation of our numerical method by solving the model presented in [41] and compare the numerical results with the results reported by [41].

Table 2. Numerical validation of $f^{″}\left(0\right)$ when $S=0$ in [41].

$*\mathit{c}/\mathit{a}.$	Present Result		Bhattacharyya et al. [41]
	First Solution	Second Solution	First Solution	Second Solution
$-0.250$	1.40224078	$-$	$1.40224051$	$-$
$-0.500$	1.49566974	$-$	$1.49566972$	$-$
$-0.625$	1.50715589	$-$	$1.50715673$	$-$
$-0.750$	1.48929822	$-$	$1.48929811$	$-$
$-1.000$	1.32881685	0	$1.32881689$	$0$
$-1.150$	1.08223113	0.11670214	$1.08223164$	$0.11670230$
$-1.200$	$0.93247330$	$0.23364972$	$0.93247277$	$0.23364910$
$-1.2465$	$0.58429940$	$0.55429554$	$0.58429146$	$0.55428565$

*c/a is the stretching/shrinking parameter in [41].

shows the comparison results, and there is a good agreement. Bhattacharyya et al. [41] employed the shooting method to solve the model, and Table 2 confirms that the bvp4c function is capable of precisely solving the boundary value problem.

Figure 2 shows the influence of the suction parameter (s) on the reduced skin friction coefficient $\left(C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}\right)$ and velocity profiles $\left(f^{\prime }\left(\eta \right)\right).$ Based on the first solution in Figure 2a, an increment in s increases the values of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}$ past a shrinking sheet. Primarily, an increment in s from 6 to 9 strengthens the impact of suction at the surface of the shrinking sheet. The act of suction encourages the laminar flow by trapping the low speed fluid molecules in the boundary layer region. This then leads to increasing of the fluid velocity past the shrinking sheet and is illuminated in Figure 2b. The increment of the fluid velocity reduces the momentum boundary layer thickness and increases the wall shear stress over the shrinking sheet. The high wall shear stress eventually increases the values of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}$ as s increases. Interestingly, the second solution in Figure 2a shows the opposite trend to the first solution, where an increment in s decreases the values of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}.$ The state of suction, which was interpreted as enhancing the fluid velocity, is now seen to decrease the fluid velocity and increase the momentum boundary layer thickness (see Figure 2). The saturated state of the shrinking sheet may be the cause of these consequences. Later, the reducing fluid velocity lowers the wall shear stress and then decreases the values of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}$ as s increases.

Figure 3 demonstrates the impact of the Deborah number (De) on $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}$ and velocity profiles. Both solutions in Figure 3a convey that an increment in De augments the values of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}$ as the sheet is shrinking. The Deborah number is used to enlighten the viscoelastic feature of a material. Here, an increment in De results in the increment of the shear thickening Sutterby fluid velocity. The valuable work of Azhar et al. [34] also reported a similar trend. The increment of the fluid velocity then enhances the wall shear stress past the shrinking sheet and increases the values of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}.$ The increment in s and De assist in delaying the flow separation significantly.

The first and second solutions in Figure 4a lead to a decrement in $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}$ when M increases from 0.5 to 1.0. The magnetic field presents in an electrically conducting fluid as an electromagnetic force in the fluid flow region, which slows the fluid moving past the shrinking sheet. This is reflected by the velocity profiles in Figure 4b, where the fluid velocity decreases when M increases. The decrement in the fluid velocity then leads to an increase of the momentum boundary layer thickness and decreases the wall shear stress past the permeable shrinking sheet. Thus, the values of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}$ decrease with the rising value of M. Unlike the shrinking case, different fluid flow behavior is perceived in the first solution when the Sutterby nanofluid flows towards a permeable stretching sheet. When the magnetic effect increases past a stretching sheet, the fluid velocity increases, although the increment is not significant. This is agreeable because the fluid flow is in the same direction as the stretching sheet and the action of the stretching sheet speeds up the fluid flow. The increment in fluid velocity reduces the momentum boundary layer thickness, increases the wall shear stress, and enhances the value of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}.$ The negative values of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}$ indicate that the stretching sheet imposes a drag force on the fluid. Moreover, the reverse flow is observed through the second solution’s presence when the permeable sheet is stretching in Figure 4a. The velocity profiles in Figure 4c support this by displaying the velocity overshoot (see the second solution profiles) when $M$ varies. Thus, it is clear that reverse flow does exist in the stretching sheet case, and this may be due to the state of the sheet where the suction intensity is weak when the effect of $M$ increases.

Figure 5 portrays the effect of the nanoparticle volume fraction or $\varphi$ on $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}$ and velocity profiles. The increment in $\varphi$ increases the values of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}$ over a permeable shrinking sheet. An increased ratio of $\varphi$ in the base fluid increases fluid viscosity, which then enhances the fluid velocity past the permeable shrinking sheet (see Figure 5b). These then affect the wall shear stress to increase and, consequently, raise the values of $C_{fx}{\mathrm{Re}}_{\overline{x}}^{1/2}.$ Velocity overshoots in the boundary layer are apparent in Figure 2b, Figure 3b, Figure 4b and Figure 5b. These velocity overshoots near the permeable shrinking sheet indicate that the fluid velocity is higher than the shrinking sheet’s velocity [42].

Table 3. The effect of the radiation parameter (Rd) on the reduced local Nusselt number when $M=0.5,m_0=1.5,\mathrm{Pr}=4.36,\varphi =0.02,s=7$, and $De=1.5$ as $\epsilon$ varies.

$\mathit{\epsilon }$	Radiation Parameter (Rd)				${\mathbf{Re}}_{\overline{\mathit{x}}}^{-1/2}\mathit{N}{\mathit{u}}_{\overline{\mathit{x}}}{\mathit{x}}^{-1/5}$
					First Solution	Second Solution
−1.5	0.5				−1846.311663	−1846.31001
−3.5					−1845.804534	−1845.803095
−5.5					−1845.296838	−1845.296378
−6.5					−1845.042742	−1845.042486
−1.5		1.2			−1846.051576	−1846.046864
−3.5					−1845.195172	−1845.191067
−5.5					−1844.337149	−1844.335834
−6.5					−1843.907426	−1843.906694
−1.5				3.5	−1845.128211	−1845.099751
−3.5					−1844.064298	−1844.051593
−5.5					−1842.553219	−1842.549137
−6.5					−1841.795465	−1841.793188

exhibits the effect of the radiation parameter (Rd) on the reduced local Nusselt number $\left({\mathrm{Re}}_{\overline{x}}^{-1/2}N{u}_{\overline{x}}{x}^{-2/5}\right)$ over the permeable shrinking surface. Both solutions allude to the enhancement of ${\mathrm{Re}}_{\overline{x}}^{-1/2}N{u}_{\overline{x}}{x}^{-2/5}$ when the impact of radiation grows in the fluid flow region. An increment in Rd hints at the release of energy in the form of heat from the fluid flow and decreases the fluid temperature profile. Thus, the thermal boundary layer becomes thinner and the wall heat flux increases. The depreciation in the thermal conductivity induces an increase in the rate of heat transfer or $\left({\mathrm{Re}}_{\overline{x}}^{-1/2}N{u}_{\overline{x}}{x}^{-2/5}\right).$ Furthermore, the magnetic parameter and the nanoparticle volume fraction have minimal effect in delaying flow separation. This is evident by the critical values $\left({\epsilon }_0\right)$, as shown in Figure 4a and Figure 5a.

The results of the stability analysis are presented in

Table 4. Lowest eigenvalues $\left({\gamma }_1\right)$ when $Rd=1.2,m_0=1.5,\mathrm{Pr}=4.36,De=1.5,s=7,M=0.5$, and $\varphi =0.02$ as $\epsilon$ varies.

ε	γ1
	First Solution	Second Solution
−6.8	0.5844	−0.4163
−6.82	0.2842	−0.1794
−6.822	0.2339	−0.1354
−6.8250	0.1128	−0.0237
−6.82520	0.0961	−0.0077
−6.825250	0.0911	−0.0028

. The first solution achieves the positive eigenvalues while the second solution attains the negative eigenvalues. Based on the signs of eigenvalues, one can say that the positive eigenvalues specify the first solution as a stable solution; the stable solution can be understood as feasible and able to overcome the growth of an initially given disturbance. Furthermore, the negative eigenvalues reveal the second solution as an unstable solution associated with flow separation. The second solution promotes the growth of an initially given disturbance and hence achieves the negative eigenvalue. However, it is vital to identify and verify the stability of non-unique solutions so that the variety of possibilities of fluid flow behavior can be predicted.

Figure 6 depicts the streamlines of the Sutterby fluid under a number of settings. In particular, Figure 6a shows the streamlines when the sheet is impermeable and stretching at the rate of 1.4, while Figure 6b illustrates the streamlines when the sheet is impermeable and shrinking. The reverse flow in Figure 6b is noticeable and proves that the shrinking sheet’s state instigates the reverse flow. Next, the streamlines for the fluid flow under the suction influence can be examined (see Figure 6c,d). The reverse flow is now absent past the permeable shrinking sheet (Figure 6d). Thus, it is proved that mass suction succeeds in sustaining the laminar boundary layer flow over a shrinking surface. Figure 6e,f shows the behavior of fluid flow when the rate of stretching or shrinking increases; the fluid pattern being pulled at the surface of the sheet is clear, and again the reverse flow is absent.

7. Conclusions

The present numerical investigation aimed to reveal the Sutterby nanofluid fluid flow and heat transfer over a permeable stretching/shrinking surface together with the effects of thermal radiation and magnetohydrodynamics (MHD). The appropriate form of the similarity transformations for the present flow problem was derived using scaling group analysis. The newly derived similarity variables then transformed the mathematical model into a more straightforward form to solve the boundary value problem utilizing the solver function bvp4c in the MATLAB software. The significant results are summarized as follows:

• An increment in the suction parameter, the Deborah number, and the nanoparticle volume fraction delay flow separation.
• The dominance of the magnetic parameter in the fluid flow regime accelerates flow separation.
• Non-unique solutions are observed when governing parameters, such as the suction parameter, the Deborah number, the magnetic number, the radiation parameter, and the nanoparticle volume fraction, vary.
• The increment in the radiation parameter slightly enhances the convective heat transfer rate past a permeable shrinking sheet.
• Stability analysis elucidates the first solution as a stable solution, and the second solution as an unstable solution.