Radiative MHD Sutterby Nanoﬂuid Flow Past a Moving Sheet: Scaling Group Analysis

: The present theoretical work endeavors to solve the Sutterby nanoﬂuid ﬂow and heat transfer problem over a permeable moving sheet, together with the presence of thermal radiation and magnetohydrodynamics (MHD). The ﬂuid ﬂow 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 di ﬀ erent numerical solutions. One of the solutions is stable and physically reliable, while the other solution is unstable and is associated with ﬂow separation. An increased e ﬀ ect of the thermal radiation improves the rate of convective heat transfer past the permeable shrinking sheet.


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 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 0 m  and 0, m  respectively. By taking the velocity field as     , , , , u x y v x y    V and under the assumptions mentioned earlier, the governing boundary layer equations in the dimensional form can be formed as follows [36]:  Sutterby [1,2] introduced the constitutive law for the Sutterby fluid by expressing the Cauchy stress tensor (T) as: where p is the pressure, I is the identity vector, and S is the extra stress tensor which can be defined as follows [21]: Here, µ 0 is the viscosity at low shear rates, E is the material time constant, . γ = tr(A 1 ) 2 /2 is the second invariant strain tensor, A 1 is the first order Rivlin-Erickson tensor or deformation rate tensor which is defined as A 1 = (∇V) + (∇V) 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 V = [u(x, y), v(x, y)], and under the assumptions mentioned earlier, the governing boundary layer equations in the dimensional form can be formed as follows [36]: where u and v denote velocity components in the x and y directions, respectively, µ n f is the dynamic viscosity of the nanofluid, σ is the electrical conductivity, B 0 is the magnetic field strength, ρ n f is the density of the nanofluid, σ 1 is the Stefan Boltzmann constant, k 1 is the Rosseland mean absorption coefficient, k n f is the thermal conductivity of the nanofluid, and C p n f 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]: where φ denotes the nanoparticle volume fraction, µ b f denotes the dynamic viscosity of the base fluid, α n f denotes the thermal diffusivity of the nanofluid, k b f 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 ρ b f and ρ 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 displays the specific values for the respective thermophysical features of n-Hexane and magnetite nanofluid [38].

Non-Dimensionalization of the Governing Equations
Considering the following the non-dimensional variables: where u 0 is the characteristic velocity and introducing the stream function ψ, which can be defined by u = ∂ψ ∂y and v = − ∂ψ ∂x , Equations (4) and (5) become: ∂ψ ∂y ∂ψ ∂y Pr Rd Pr with the corresponding boundary conditions: while satisfying the continuity equation of Equation (3). In Equations (9) and (11), M = is the Deborah number, φ is the nanoparticle volume fraction, and terms A 1 , A 2 , A 3 , and A 4 are expressed as: The functions u w , v w and T w (x) are assumed to be in the following form to ensure that similarity solution exists: where u 1 is the reference velocity, v 1 is the normal reference velocity, and T 0 is the reference temperature.

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: where ω i are constants to be determined in which i = 1, . . . 8. The transformation G is the transformation point which transforms the (x, y, ψ, σ, θ, u e , u 1 , m,) coordinates to the new coordinates x * , y * , ψ * , σ * , θ * , u * e , u * 1 , m * . Next, the substitution of (14) into Equations (9)-(11) yields the following expressions: To retain the invariance of the system under G, the parameters defined in Equation (14), the following relations must hold: From the boundary conditions of Equations (17), we also obtain the following relations among the parameters: The absolute invariant can be determined by eliminating the parameter G of the group and hence Equations (18) and (19) provide the following expressions: From Equations (13), (14), and (20), we achieve the absolute invariants under the group G similarity transformations as follows: 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): with the associated boundary conditions: Here ε = (u 1 ) 0 /u 0 is the stretching/shrinking parameter, where ε > 0 indicates the stretching sheet, ε = 0 specifies the stationary sheet, and ε < 0 represents the state of shrinking sheet. Furthermore, 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 (u e ) 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 C f x and the local Nusselt number (Nu x ) which are defined as follows: where τ 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]: The reduced skin friction coefficient Re 1/2 x C f x x 1/5 and the local Nusselt number Re −1/2 x Nu x x −1/5 can be obtained using the similarity transformations of Equation (21) and the expressions in Equations (25) and (26) as follows: where Re x = u e x/ν b f denotes the local Reynolds number.

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: ∂T ∂t with the boundary conditions of Equation (6). Then, we introduce the dimensional time variable, t with a new similarity variable τ = τ 0 /x 4/5 through scaling group analysis. The new similarity transformation is given as: Employing (31) in the dimensionless form of Equations (28)-(30) and (6) gives the following system of equations: with the boundary conditions: It is assumed that the solutions of (32)- (34) are expressed by the formulas of Equation (35): where f (η) = f 0 (η) and θ(η) = θ 0 (η) 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 γ, and F(η, τ) and G(η, τ) are relatively small compared to the steady state solutions ( f 0 (η) and θ 0 (η)). The substitution of Equation (35) into Equations (32)-(34) gives the following system: subject to the boundary conditions: Referring to Merkin [39,40], τ → 0 is fixed to examine the stability of the steady state boundary layer flow. Thus, F = F 0 (η) and G = G 0 (η) in Equations (37)-(39), yielding the following linearized eigenvalue problem: with the boundary conditions: 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 (∞) = 0 is substituted with F 0 (0) = 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.

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 × 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 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 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.  x . 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 f x Re 1/2 x as s increases. Re fx x C 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 12 Re fx x C as s increases.
(a) (b)   Re fx x C Velocity overshoots in the boundary layer are apparent in Figures 2b, 3b, 4b, and 5b. These velocity overshoots near the permeable shrinking sheet indicate that the fluid velocity is higher than the shrinking sheet's velocity [42].  x over a permeable shrinking sheet. An increased ratio of φ 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 f x Re 1/2 x . Velocity overshoots in the boundary layer are apparent in Figures 2b, 3b, 4b and 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 results of the stability analysis are presented in Table 4. 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 results of the stability analysis are presented in Table 4. 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.   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.

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

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.