1 Introduction

The non-Newtonian fluid has intrigued many scholars due to its enormous significance in nature and engineering systems, for instance, electronic chips, polymer depolarization, boiling, fermentation, composite processing, and bioengineering. In 1867, Maxwell^[1] proposed a new kind of fluid, known as the Maxwell fluid. It is well-known that this type of fluid is the most basic subclass in the non-Newtonian fluid. The Maxwell fluid model estimates effects on relaxation time and is capable of estimating the shear-thinning liquid. This model also removes the complicated effects of shear-dependent viscosity, allowing researchers to focus solely on the role of fluid elasticity on boundary layer parameters^[2]. Wang and Tan^[3] investigated the stability performance of the Maxwell fluid model, and observed that the relaxation time boosted the system's instability.

The exploration of the Maxwell fluid model was broadened by Nadeem et al.^[4] who incorporated the effects of nanoparticles in the non-Newtonian Maxwell fluid past a stretching sheet. The results showed that the inclusion of nanoparticles revealed a significant effect on the heat transfer and behaviors of the working fluid. The experimental work of Binetti et al.^[5] proved that the sodium alginate (SA) exhibits the linear viscoelastic feature, which is an example of the Maxwell fluid. Next, Pawar and Sunnapwar^[6] undertook an investigation on the heat transfer performance of non-Newtonian fluids (including the SA) in a helically coiled heat exchanger. Their study has provided the thermophysical properties data of the SA comprehensively. The SA has also been attempted as nanoparticles^[7] and as base fluids^[8], which benefits the applications in agriculture sectors.

In recent years, growing attention has been paid to classifying the nanoparticles' properties in a Newtonian hybrid nanofluid. However, limited studies have been reported on non-Newtonian hybrid nanofluid properties. Bahiraei et al.^[9] performed an experimental study on the hydraulic and thermal features in a mini channel heat exchanger. The obtained results reflect a favorable outlook for using the non-Newtonian hybrid nanofluid in mini heat exchangers. Al-Rashed et al.^[10] presented a numerical analysis of non-Newtonian hybrid nanofluid considering the thermodynamic properties in a microchannel heat sink. The results demonstrate that by increasing the concentration of nanoparticles, the thermal resistance decreases. Other studies in which non-Newtonian fluid was used to assess thermal efficiency can be reviewed here^[11-13]. It is also worth mentioning here some good papers that discussed the respective problems^[14-16].

Radiation parameters are one of the essential control factors for heat and fluid flow in the thermal system involving high temperatures. In particular, the effects of thermal radiation on the development of stable equipment, gas turbines, nuclear plants, satellites, missiles, and several complex conversion systems are very important^[17-18]. Accordingly, Cess^[19] was the first researcher who studied the heat transfer interaction between thermal reaction and free convection by employing the perturbation technique. The remarkable work has been widened by Arpaci^[20] in a heated vertical plate towards a stagnant radiating gas. Since then, countless research works have been carried out on concerns with thermal radiation, such as Agbaje et al.^[21], who performed a numerical analysis on the unsteady flow of Powell-Eyring nanofluid in the presence of thermal radiation and heat generation. The results revealed that the heat generation and thermal radiation parameters increased the thermal boundary layer thickness and fluid temperature.

In a nutshell, non-Newtonian fluid dynamics has gained much importance in recent decades due to its numerous practical applications. The typical fluid flow characteristics in the industry, such as polymeric liquids, biological fluids, and motor oils, can also be defined via the non-Newtonian fluid model. A review of the previously listed literature, however, addresses a variety of limitations and shortcomings. This research work aims to construct a mathematical model driven by the initial analysis to assess the numerical solutions of the unsteady flow in the Maxwell hybrid nanofluids with thermal radiation impact. The thermophysical characteristics of hybrid nanofluids are adopted from Takabi and Salehi^[22] and Ghalambaz et al.^[23]. It is important to note that an effective numerical procedure is crucial in solving this mathematical model. Hence, in this study, the bvp4c scheme in the MATLAB software is prompted to solve the resulting ordinary differential equations. This solver is built by utilizing a collocation approach with the fourth-order precision. By employing a pair of arbitrary initial guesses, the bvp4c solver can efficiently anticipate solutions. Also, the method described above has successfully provided multiple solutions in the current study. Moreover, previous literature demonstrated that the bvp4c solver is accurate in addressing boundary layer problems, as it is able to withstand and compete with other numerical methods such as the Keller-box method and the shooting method^[24-26].

2 Mathematical model

In this study, we deliberate the unsteady upper convected Maxwell hybrid nanofluid over a permeable stretching/shrinking sheet and heat transfer in two-dimensional laminar boundary layer flow, which starts at time t=0. The sheet moves with the non-uniform velocity u_w (x, t), where (x, y) are the Cartesian coordinates, and the x- and y-axes are assigned along the surface and normal to it while the flow is at y≥ 0 (see Fig. 1). Also, the surface temperature T_w and ambient temperature T_∞ are assumed to be constant. All these assumptions apply to geophysical problems such as double diffusive convection or salt finger evolution, reactor design in chemical processing, and turbidity currents^[27-28]. Therefore, the governing equations of the respecting assumptions can be defined as^{[18, 29]}

(1)

(2)

(3)

subject to

(4)

Here, u is the velocity along the x-direction and v in the y-direction, v_w (t) represents the mass flux velocity with v_w (t) < 0 for suction and v_w (t)>0 for injection, respectively, and α denotes the constant having dimension time t^-1 such as (αt < 1, αt≥ 0), a denotes a positive constant, q_r is the radiation, μ_hnf is the Al₂O₃-Cu/SA dynamic viscosity, k_hnf is the Al₂O₃-Cu/SA heat/thermal conductivity, ρ_hnf is the Al₂O₃-Cu/SA density, T is the temperature of Al₂O₃-Cu/SA, and (ρc_p)_hnf is the Al₂O₃-Cu/SA heat capacity. Table 1 shows the hybrid nanofluid thermophysical properties as defined by Takabi and Salehi^[22] and Ghalambaz et al.^[23]. Here, ϕ is the volume fraction of nanoparticles, ρ_f denotes the SA density, ρ_s is the hybrid nanoparticle density, k_f and k_s indicate the thermal conductivity of SA and hybrid nanoparticles, respectively, and finally c_p is the heat capacity pressure. At the same time, the thermophysical properties of Cu, Al₂O₃ nanoparticles and also base fluid, SA, are demonstrated in Table 2 (see Ref. [30]).

The radiative heat flux q_r is presented by the Rosseland^[31] approximation,

(5)

where k^* is the coefficient of mean absorption, and σ^* is the Stefan-Boltzmann constant. We presume that within the flow, the temperature difference T⁴ can be extended in Taylor's series. Hence, enlarging T⁴ about T_∞ while disregarding the terms of the higher-order, we obtain

(6)

Thus,

(7)

Substituting Eq. (7) into Eq. (3), we get

(8)

By considering the following dimensionless variables^[29]:

(9)

with the stream function ψ, which is defined as and , Eq. (1) is fully satisfied. Then, we have

(10)

so that

(11)

In the above equation, S is a parameter of the mass flux, where S=0 denotes an impermeable surface, and S < 0 and S>0 correspond to injection and suction, respectively. Substituting (9) into Eqs. (2) and (8), we get the following ordinary (similarity) differential equations:

(12)

(13)

subject to

(14)

where primes denote differentiation with respect to η, Pr=ν_f/α_f is the Prandtl number, λ is the constant stretching/shrinking parameter with λ >0 for the stretching sheet, λ < 0 for the shrinking sheet, and λ =0 for the static sheet. The unsteadiness parameter is defined by β =a/α with β < 0 indicating the decelerating flows, and β >0 indicating the accelerating flow, while β =0 indicates the steady-state flow, and K=ak₀ is the Maxwell parameter. Now, we describe C_f (skin friction coefficient) and Nu_x (local Nusselt number) as follows:

(15)

We note that τ_w and q_w are the shear stress and the heat flux from the plate, respectively. Thus, we have

(16)

Utilizing Eqs. (9) and (16), we get

(17)

where Re_x=u_w(x, t)x/ν_f is the local Reynolds number, and R_d=4σ^* T_∞³/(k_fk^*) is the radiation parameter.

3 Stability analysis

Since several solutions in the boundary value problem (12)-(13) have been demonstrated, an examination of solution stability is carried out. This technique has been frequently adopted by a number of researchers due to its validity^[32-34]. By perceiving the efforts of Merkin^[35] and Merrill et al.^[36], we instigate a new transformation variable τ under the unsteady-state query in the following way^[37-38]:

(18)

Substituting Eq. (15) into Eqs. (12)-(13), the following equations are secured:

(19)

(20)

with respect to

(21)

To scrutinize the steady flow solutions stability, where f(η)=f₀ (η) and θ (η)=θ₀ (η), we describe f(η, τ)=f₀(η)+e^-ωτF(η) and θ (η, τ)=θ₀ (η)+e^-ωτH(η), which is developed by Weidman et al.^[39]. Here, F(η) and H(η) are relatively small to f₀ (η) and θ₀ (η), while ω is the undetermined eigenvalue parameter. The eigenvalue problem (19)-(20) leads to an infinite number of eigenvalues ω₁ < ω₂ < ω₃⋯ that disclose a steady flow and initial deterioration when ω₁ is positive. On the other hand, once ω₁ is negative, an initial development of perturbations will be discovered. This finding exposes the occurrence of unstable flow. By substituting Eq. (21) into Eqs. (19)-(21), the equations below are generated,

(22)

(23)

(24)

Now, we implement the steady-state flow solutions f₀(η) and θ₀ (η) via τ→0. As a result, the linearized eigenvalue problem is determined as follows:

(25)

(26)

(27)

The possible eigenvalues could be deliberated via relaxing a boundary condition^[40], where the value of F'(η)→0 is fixed as η→∞, and thus the linearized eigenvalue problems (25)-(26) are exposed as F"(0)=1 when ω₁ is set up.

4 Interpretation of results

The resulting ordinary differential equations presented in Eqs. (12)-(14) are clarified in the MATLAB software^[41] by utilizing the bvp4c method. In this technique, the essential approximation of boundary layer thickness is critical while discovering multiple solutions. Before achieving the required outcome, several attempts need to be performed in offering a significant initial prediction. Further, the reliability of the findings is measured with Sharidan et al.^[42], who employed the Keller-box method and Chamka et al.^[43], together with Madhu et al.^[29], who utilized the finite-difference method. By establishing an explicit agreement with previous literature, the significance of the existing numerical approach, i.e., the bvp4c function, is undeniable as accessible in Table 3. Thus, we suggest that the proposed conceptual model thoroughly examine heat transfer and fluid flow activity with considerable certainty.

Figures 2-15 unveil the existence of the dual solution for preferred proportions of the controlling parameters when the stretching/shrinking parameter (λ) is varied. In the entire study, the dual solutions are seen as the sheet shrinks (λ < 0) within a certain limit of a critical point, i.e., λ_c. When λ < λ_c, no solution is noticed as it expresses the separation of the boundary layer. It is also worth noting here that the counter-argument of which solutions ultimately take place in a real-life application depends on flow stability, which demands comprehensive examination. Furthermore, certain limiting parameters are utilized to assure the accuracy of the outputs, and they are fixed to the extents listed here; 0.0≤ ϕ₂ ≤ 0.3, 2.0≤ S≤ 2.4, 0.03≤ K≤ 0.08, -0.5≤ β ≤ -1.5, and 1≤ R_d≤ 3. During the analytical execution of this investigation, the Prandtl number is assigned to 6.5, omitting any similarities with previous arguments.

Figures 2 and 3 expose the effects of assorted values in the suction parameter S toward λ past a stretching/shrinking surface. The characteristic of f"(0) in Al₂O₃-Cu/SA is described in Fig. 2, and it indicates that, in the first solution, an improvement in S would substantially upsurge f"(0). By comparison, another solution reveals a declining tendency of f"(0) with the addition of S. The occurrence of suction on a stretching/shrinking sheet potentially promotes the stabilization of the boundary layer. Besides, the suction parameter lowers the body pressure in such a decelerating flow, thereby augmenting the velocity gradient and reducing the thickness of the stretching/shrinking sheet boundary layer by evacuating the fluid with low momentum alongside the surface. This is proven in Fig. 4, which displays the velocity profile f'(η) with multiple values of S. Meanwhile, the first and second solutions provide an enhancement in -θ'(0) as S escalates through the Al₂O₃-Cu/SA stretching/shrinking surface, which is clarified in Fig. 3. It is worth noting that the suction action allows Maxwell hybrid nanofluid particles to occupy the surface as it stretches or shrinks and boost up the rate of heat transfer physically. The temperature distribution profile θ (η) confirms this argument by showing a downward trend in both solutions, where suction is imposed in the stretching/shrinking sheet, as displayed in Fig. 5. In general, the momentum and temperature distribution profiles f'(η) and θ (η) satisfy the far-field boundary conditions (14) whenever η_∞ =12 is executed.

Figures 6 and 7 inaugurate the effects of Maxwell parameter K on f"(0) and -θ'(0) in the Maxwell hybrid nanofluid, respectively. The consequence of rising K values is to reduce the velocity field and thereby decrease the thickness of the boundary layer, as presented in Fig. 8. From a physical perspective, fluid tends to rest when the shear stress is eliminated. This phenomenon can be observed in a variety of polymeric liquids that defy the Newtonian fluid model. Also, a high value of the Maxwell parameter may create a retarding force in the flow between two adjacent layers. This leads to a decline in the velocity and the thickness of the boundary layer, which can be seen in Fig. 8. Figure 7 displays a diminishing response of -θ'(0) when K varies in the first solution. Concisely, the rate of heat transfer decreases as K improves in a Maxwell hybrid nanofluid. This observation suggests that the addition of the Maxwell parameter in the boundary layer may minimize the heat transfer rate of the thermal boundary layer, thus increasing the fluid's temperature. This is in accord with the temperature profile trend at the stretching/shrinking surface of the Maxwell hybrid nanofluid Al₂O₃-Cu/SA as disclosed in Fig. 9. Figure 9 confirms that by extending the Maxwell parameter value, the temperature profile increases since the thickening of the thermal boundary layer enhances the elasticity stress parameter. Additionally, improving the value of K in the Maxwell hybrid nanofluid is proven to boost up characteristics of the first solution; however, diminishes in other solution of the temperature distribution profile.

The effects of the unsteadiness parameter β toward λ when β shifts from -0.5 to -1.5 are demonstrated in Figs. 10-11. The Maxwell hybrid nanofluid Al₂O₃-Cu/SA characteristic is depicted in fig. 10 with respect to f"(0) when β is reduced. the reduction of β leads to the expansion of the boundary layer thickness and subsequently declines the velocity gradient, and hence f"(0) diminishes. the presence of the nanoparticle volume fraction could also initiate the decline of f"(0) because of the rise in the maxwell hybrid nanofluid viscosity. moreover, according to the generated results in fig. 11, -θ'(0) is increased in both solutions as β decreases. based on the present results, the authors can deduce that a reduction value of the unsteadiness parameter promotes the rise of the heat transfer rate significantly. even so, if multiple control parameters are taken into account, the authors would also like to claim that those effects may vary.

Figures 12-13 present the responses of f"(0) and -θ'(0) in the Maxwell hybrid nanofluid when ϕ₂=0.01, 0.02, and 0.03. Figure 12 expresses that an augmentation in the nanoparticle volume fraction ϕ₂ (Cu) is proven to intensify the f"(0) values. The intensification of ϕ₂ improves the Maxwell hybrid nanofluid viscosity, which successively boosts the fluid velocity. While f"(0) increases, the frictional drag exerted increases, potentially slowing the separation of the boundary layer as the surface shrinks. In contrast, the downward trend of f"(0) is revealed when ϕ₂ improves, and the value of reduced skin friction coefficient progressively becomes negative. The negative value of f"(0) elucidates that the Maxwell hybrid nanofluid imposes the drag force away from the sheet. Figure 12 also emphasizes that when λ=0 (static sheet), f"(0)=0, which reflects a lack of frictional drag on the surface of the static sheet. Figure 13 depicts the behavior of declining thermal performance or -θ'(0) towards a shrinking sheet, and this behavior can be seen in every solution. In other words, as ϕ₂ rises, the heat transmission performance decreases. Theoretically, as the mass of the Maxwell hybrid nanoparticles increases, the temperature of the sheet rises; as a result, the efficiency of thermal conductivity drops, affecting the performance of heat transfer. This could be due to radiation parameter effects that raising the temperature across the surface and disrupting nanoparticle interactions, result in lowering thermal conductivities.

The effects of the radiation parameter R_d in the Maxwell hybrid nanofluid Al₂O₃-Cu/SA are shown in Figs. 14 and 15. Figure 14 presents the variants of -θ'(0) for several values of R_d where R_d=1.0, 2.0, and 3.0, while Fig. 15 displays the temperature profile distribution θ (η). The thermal radiation parameter is discovered to increase the temperature profile distribution. Substantially, -θ' (0) boosts when the thermal radiation parameter is intensified attributing to the conversion of thermal energy. This is due to the development and transmission of more heat into the flow, which helps increase the momentum boundary layer thickness. The same findings were found in the preceding literature^[18]. On another note, Fig. 15 displays the increment of θ (η) when R_d improves to the extent that the energy conversion rate over the non-Newtonian Maxwell hybrid nanofluid is enlarged. The observations in the previous analysis (see Ref. [44]) indicated the same effects as the current findings. Up to now, it has been demonstrated that the influence of thermal radiation can raise the temperature of a Maxwell hybrid nanofluid in a stretching/shrinking surface prior to the listed condition.

A stability analysis was further conducted by utilizing the bvp4c application in the MATLAB software. The smallest eigenvalues ω₁ for certain values of λ when ϕ₁=ϕ₂=0.01, S=2.2, K=0.05, β =1.0, R_d=1.0, and Pr=6.5 are given in Table 4. The first solution is usually denoted as reliable because it fulfils the far-field boundary condition. By conducting a stability analysis, we are able to confirm the accurate and reliable solutions convincingly. In the stability analysis procedure, the smallest eigenvalue ω₁ indicates the characteristics of the dual solutions. When ε₁ is positive, the flow is considered stable since the solutions comply with the stabilizing property where an initial decay is permitted. Meanwhile, when ω₁ appears negative, the flow is represented unreliable as it evokes an initial extension of disturbances. The first solution is found to be reliable, while the alternative solution is unreliable, as witnessed in Table 4.

5 Conclusions

A mathematical analysis of the unsteady Maxwell hybrid nanofluid over a stretching/shrinking surface with the impact of thermal radiation is validated. The Maxwell hybrid nanofluid is chosen as the operating fluid since this model is proposed to describe the behavior of viscoelastic fluids and succeed in modeling the polymeric liquids. The present work is an example of a numerical approach to a thermal extrusion manufacturing process, such as drawing plastic films or extrusion of a polymer sheet from a dye. The findings indicate that the existence of dual solutions across the Maxwell hybrid nanofluid Al₂O₃-Cu/SA flow is factual in specific controlling parameters. The suction action permits Maxwell hybrid nanofluid molecules to occupy the sheet, resulting in the increment of the heat transfer rate, physically. When the thermal radiation parameter is increased, the heat transfer rate boosts due to the presence of thermal energy conversion. The augmentation of the nanoparticle volume fraction raises the Maxwell hybrid nanofluid viscosity, resulting in an upsurge in fluid velocity. A reduction value of the unsteadiness parameter that denotes the decelerating flow promotes the heat transfer rate significantly. Meanwhile, an increment in the Maxwell parameter value causes a decrement in the heat transfer rate due to the retarding force in the flow. Finally, the stability analysis guarantees the first solution's consistency, while the second solution is declared unstable.