Taguchi–grey relational analysis and response surface methodology for heat transfer optimization in tri-hybrid nanofluid flow with slip and heat source/sink effects

1 Introduction

As the interest on nanofluids, hybrid nanofluids and THNFs as advanced heat-transfer media is increasing, the application of systematic optimization techniques has become more important to explore the effect of multiple operating parameters on thermal performance. RSM has been recognized as an effective tool for process modeling, optimization and performance evaluation among the available experimental and statistical approaches. RSM is a set of mathematical and statistical methods to define the functional relationship between a number of independent variables and one or more response variables. The main aim of RSM is to find the optimum operating conditions and improve the desired responses [1–3]. The foundation for RSM was laid by Box et al. [4], who presented the methodology as a framework for experimental design, empirical model building, and the study of interactions among process variables. Their work showed the potential of RSM to identify operating regions that give good results with a minimum of experimental effort. Based on these principles, Box et al. [5] made important contributions, such as steepest ascent approach and Central Composite Design (CCD), which allowed the efficient estimation of second-order response surfaces and improved optimization accuracy. Further progress was reported by Eide et al. [6] who extended the use of RSM to mixture-design problems where the system responses depend upon the relative proportions of the constituent components rather than upon independent process variables. This extension has significantly extended the scope of application of RSM in fields such as chemical engineering, pharmaceutical sciences and materials technology. In recent years, RSM has been combined with modern computational and artificial intelligence methods such as genetic algorithm (GA), artificial neural network (ANN), particle swarm optimization (PSO), and machine-learning based approaches to enhance both predictive performance and optimization efficiency for complex engineering systems [7–10]. The ever-evolution of RSM has greatly expanded its usage in various fields like nanofluid engineering, optimization of thermal systems, manufacturing technologies, biotechnology, and advanced materials science. RSM is one of the most efficient and commonly used optimization techniques in scientific and engineering studies because it can consider the combined influence of several parameters and find suitable operating conditions [11, 12]. ANOVA is one of the most important steps in RSM. It is used to test the statistical significance of the model terms and to estimate the contribution of individual factors and their interactions to the observed responses. Several researchers have successfully applied RSM and ANOVA in heat-transfer and fluid-flow investigations. For example, Gupta et al. [13] used these techniques to explore the thermal performance of heat pipes working with nanofluids and shown the efficiency of these techniques in identifying influential process parameters. Pattnaik et al. [14] optimized the thermal characteristics of a radiative Carreau THNF system using RSM and developed predictive relationships among the governing variables and response functions. Raza et al. [15] also used RSM-based optimization and observed that the maximum heat-transfer effectiveness occurred at the higher concentration of nanoparticle and higher Reynolds numbers. Moreover, RSM and statistical analysis have been used in many studies to inspect complex fluid-flow and heat-transfer phenomena, which prove the versatility and reliability of these approaches in engineering applications (Jebali et al. [16], Mishra et al; [17], Wantao et al; [18], Abu Bakar et al; [19], Parida et al; [20]).

A highly efficient approach for enhancing heat transfer characteristics is without sacrificing cost-effectiveness is the application of nanofluids as heat-conducting coolants. In the recent years, there is a great interest in improving the heat transfer characteristics using THNF in different thermal engineering systems. Typically, a nanofluid is a suspension of a particular type of nanoparticles in a conventional base fluid. In this context, HNFs are obtained by dispersing two different types of nanoparticles in a base fluid, whereas THNFs are produced by adding three different types of nanoparticles in the same host fluid. The high thermal conductivity of these nanostructured materials is an important factor for the improvement of heat transfer processes. The uniform dispersion of the nanomaterials in the base fluid improves the overall thermophysical properties of the suspension and hence the thermal transport efficiency. Due to such promising properties, THNFs have attracted considerable research interest and are being progressively explored for their potential use in advanced cooling technologies, energy systems, manufacturing processes and other engineering applications. Ali et al. [21] observed that the improvement in energy transfer was more pronounced for THNFs compared to mono and HNFs in the heat sink conditions. Dey et al. [22] observed a significant enhancement of the heat-transfer coefficient and concluded that the synergistic interaction between three constituent nanoparticles allows THNFs to perform better than conventional nanofluids and HNFs, both in terms of velocity and temperature distributions. Correspondingly, Saadeh et al. [23] observed that THNF could enhance the temperature profile by 31.6% and velocity field by 27.4% in comparison with the corresponding base nanofluid. Moreover, Yasir et al. [24] stated that the enhancement of the thermal transport properties is significantly improved with the addition of ternary hybrid nanoparticles as indicated by the higher values of N u x. These results collectively highlight the potential of THNFs as highly effective heat-transfer fluids. Some more studies on the different aspects of THNF flow and heat-transfer behavior are available in the literature [25–28].

Many industrial and engineering systems such as nuclear reactors, geothermal energy extraction, thermal energy storage units, electronic cooling devices, chemical reactors and advanced technologies of thermal management have internal heat generation and absorption mechanisms. Thus, the inclusion of heat sink/source impacts in fluid-flow models has become an important research topic for the enhancement of thermal performance and control of heat-transfer processes. The practical importance of the combination of heat generation/absorption with thermal slip, stretching/shrinking surfaces, and suction conditions in the context of THNFs has attracted considerable attention in the context of modern heat-transfer applications. Some researchers have studied the effects of heat source/sink on nanofluid and THNF flows. Khan et al. [29] studied the effects of heat generation and absorption on the THNF boundary-layer flow over a stretching surface and showed its significant effect on thermal characteristics. Mandal et al. [30] investigated a Casson-Williamson THNF with heat source/sink effects and reported an improvement in the heat transfer rate. In the same way, Alam et al. [31] considered a radiative THNF under the influence of an exponential heat source/sink and revealed its substantial effect on the temperature distribution in the boundary layer (BL). Ho et al. [32] studied the heat-transfer characteristics of a mini-channel heat sink/source system using an Al₂O₃–water nanofluid and reported a considerable improvement in the average N u x at both the heating and cooling sections. Irfan et al. [33] investigated the combined effects of magnetic field and heat source/sink on Maxwell nanofluid flow and found significant variations in thermal and flow features. Further Khan et al. [34] investigated the influence of nanoparticle mass flux theory and heat generation/absorption on radiative Maxwell nanofluid flow over a stretching cylinder. The results show that the temperature profile decreases with the increase of the Maxwell fluid parameter, the Deborah number. These investigations collectively demonstrate the crucial role of heat source/sink effects on thermal transport processes and motivate future studies with more sophisticated THNF configurations.

The governing equations are numerically solved with the MATLAB bvp4c solver, however the numerical solution will only provide predictions of flow and heat transfer characteristics for specific values of the governing parameters. In the applications of practical thermal engineering, however, the designer is interested in finding the combination of operating parameters which will result in heat-transfer as high as possible, while maintaining acceptable hydrodynamic characteristics. The numerical search procedure over all the possible combinations of parameters is also very expensive since a single combination needs to be solved numerically for the boundary values. To obtain the predictive relationships between the governing parameters and the engineering responses, statistical optimization technique is used. The Taguchi method offers an efficient numerical experiments design, Grey Relational Analysis (GRA) transforms several performance measures into one optimization criterion and Response Surface Methodology (RSM) builds predictive models, quantifies the interactions of the parameters, and determines the optimal operating conditions. The statistical optimization is thus complementary to the deterministic numerical simulations, and is used to translate numerical data into practical design guidelines for engineering applications.

Based on the above studies and research gap, the present study is planned to investigate the combined effects of suction, stretching/shrinking, thermal slip, heat source/sink and trihybrid nanoparticle concentration on the flow and heat transfer characteristics of THNF. Special emphasis is given to the existence of multiple solutions and hence all possible numerical solutions are generated and then investigated using the stability analysis to identify the physically realizable stable solution branches. The nonlinear governing equations are transformed into a system of ordinary differential equations and solved numerically using the MATLAB bvp4c solver. Furthermore, a statistical analysis based on the RSM is carried out to quantify the effect of the governing physical parameters on the engineering quantities of interest. The experimental plan for ANOVA employs the Central Composite Design (CCD) for the construction of predictive models for C f x and N u x. To the best of authors' knowledge, the combined effects of suction, stretching/shrinking, thermal slip, heat source/sink and trihybrid nanoparticle concentration have not been considered together in a single framework while integrating stability analysis with RSM-based statistical modeling in any previous study. Thus, the current study presents a new and complete approach to understand and improve the heat transfer performance in THNF flows.

2 Mathematical formulation

In this work, the two-dimensional steady boundary-layer motion of THNF over a permeable shrinking/stretching sheet subjected to a power-law heat flux condition is examined. The formulation further accounts for thermal slip, wall suction, and volumetric heat generation/absorption effects. A schematic depiction of the flow relationship together with the coordinate framework is presented in Figure 1. The governing mathematical model is established under the subsequent assumptions:

• The carrier fluid and dispersed nanoparticles are considered to be in thermal equilibrium, and the THNF mixture is treated as homogeneous.
• The surface undergoes stretching or shrinking with the velocity u w x = α x λ, where α denotes a positive constant and λ is the stretching/shrinking parameter. A shrinking sheet is represented by (λ < 0), whereas (λ > 0) characterizes a stretching sheet.
• The sheet is permeable and is subjected to a specified mass transfer velocity v w.
• A heat source/sink parameter is included to account for the generation or absorption of thermal energy within the fluid domain.
• The wall is subjected to a power-law heat flux of the form q w = q 0 x 2 where q 0 is the constant heat flux coefficient, and 2 is the heat flux power-law index.

The governing continuity, momentum, and energy equations under the above assumptions are written as follows [35, 36]: along with boundary conditions (BCs)

v y = − u x (1)

u u x ‐ μ thnf ρ thnf u yy + v u y = 0 (2)

u T x ‐ k thnf ρ C p thnf T yy + v T y ‐ Q ρ C p thnf T ‐ T ∞ = 0 (3)

v − v w = 0 , u − u w x = 0 , T w + T S l i p − T = 0 a t y = 0 T → T ∞ , u → 0 , a s y → ∞ (4)

The temperature of THNF is denoted by T, and the free-stream temperature is represented by T ∞. The thermal slip at the boundary is given by T Slip = D 1 ∂ T ∂ y which is a first-order slip condition, where D 1 is the thermal slip coefficient. Q is the heat absorption/generation parameter. The base medium is assumed to be water (H 2 O) as the working fluid, with the presence of alumina (A l 2 O 3), copper (Cu), and copper oxide (CuO) nanoparticles. Thermophysical parameters with their relationship with the fluid system are shown in Figure 2, and the thermophysical properties of A l 2 O 3 (ϕ 1), Cu (ϕ 2), and CuO (ϕ 3) are tabulated in Table 1. Following the formulation of Lund et al. [37], we introduce the following similarity transformations:

u = α x f ′ η , v = − α ν f 1 2 f η , η = y α / ν f 1 2 , θ η = T − T ∞ T w x − T ∞ (5)

TABLE 1. Intrinsic thermophysical properties of fluid and nanoparticles [42].

The surface mass transfer is defined through v w = − α ν f 1 2 S, where (S) is a dimensionless mass flux parameter. A positive value of (S) indicates fluid withdrawal (suction), whereas a negative value corresponds to fluid blowing (injection). The wall temperature condition is T w = q w + T ∞ . By using the following values, which are derived from the help of Equation 5 u x = α f ′; u u x = α 2 x f ′ 2; u y = α x ν f f ″; v u y = − α 2 x f f ″; u y y = ν f α 2 x f ‴; T x = 2 q 0 x θ; u T x = 2 α q 0 x 2 f ′ θ; T y = q 0 x 2 ν f α θ ′; v T y = − α q 0 x 2 f θ ′; and T y y = ν f α q 0 x 2 θ ″. Now, we substitute the values and Equation 5 into Equations 2 and 3 with BCs Equation 4, we attain: subject to BCs:

f ‴ + ρ thnf / ρ f μ thnf / μ f f f ″ ‐ f ′ 2 = 0 (6)

k thnf / k f Pr ρ c p thnf / ρ c p f θ ″ + f θ ′ ‐ 2 f ′ θ + β ρ c p thnf / ρ c p f θ = 0 (7)

f 0 = S , f ′ 0 = λ , θ 0 = 1 + δ T θ ′ 0 f ′ η → 0 , θ η → 0 as η → ∞ (8)

The prime in the above equations denotes the derivative with respect to η. The parameter λ is the rate at which the surface is stretched or shrunk. The heat generation/absorption parameter is called β and is defined as β = Q α ρ C p f where β > 0 corresponds to internal heat generation (heat source), and β < 0 corresponds to heat absorption (heat sink). Also, δ t = D 2 α ϑ f is the thermal slip parameter. The C f x and N u x are expressed as:

C f x = μ t h n f ρ f u w 2 ∂ u ∂ y y = 0 , N u x = − x k t h n f k f T w x − T ∞ ∂ T ∂ y y = 0 (9)

Using Equation 5 in Equation 9, we get:where R e x = u w x x ν f is the local Reynolds number.

Re x 1 / 2 C f = μ thnf μ f f ″ 0 , Re x ‐ 1 / 2 N u x = ‐ k thnf / k f θ ′ 0 (10)

3 Stability analysis

The boundary value problem defined by Equations 6 and 7 has multiple solution branches, and therefore their stability must be studied. Following the approach suggested by Merkin [38], the transformations in Equation 4 are rewritten by introducing a new variable, τ, as:

u = α x ∂ f ∂ η η , τ , v = ‐ α ν f 1 2 f ( η , τ ) , η = y α ν f 1 2 θ η , τ = T ‐ T ∞ T w x ‐ T ∞ , τ = α t (11)

Consequently, Equations 6 and 7 serve as the governing framework for

μ thnf / μ f ρ thnf / ρ f ∂ 3 f ∂ η 3 + f ∂ 2 f ∂ η 2 ‐ ∂ f ∂ η 2 ‐ ∂ 2 f ∂ η ∂ τ = 0 , (12)

k thnf / k f Pr ρ c p thnf / ρ c p f ∂ 2 θ ∂ η 2 + f ∂ θ ∂ η ‐ 2 ∂ f ∂ η θ + β ρ c p thnf / ρ c p f θ ‐ ∂ θ ∂ τ = 0 , (13)

Furthermore, Equation 8 yields the following set of conditions:

λ − ∂ f 0 , τ ∂ η = 0 , S − f 0 , τ = 0 , θ 0 , τ = 1 + δ T ∂ θ 0 , τ ∂ η , ∂ f η , τ ∂ η → 0 , θ η , τ → 0 as η → ∞ (14)

Subsequently, following the methodology presented by Weidman et al. [39], the following expression is obtained:

f η , τ − f 0 η = F η , τ e − ξ τ , θ η , τ − θ 0 η = I η , τ e − ξ τ (15)

In the stability analysis, F η and I η denote the disturbance modes corresponding to f 0 η and θ 0 η, respectively, and (ξ) represents the associated eigenvalue. Substitution of Equation 15 into Equations 12 and 13, followed by linearization of the resulting equations, yields the following system: subject to BCs:

F ‴ + ρ thnf / ρ f μ thnf / μ f f 0 ″ F + f 0 F ″ ‐ 2 f 0 ′ ‐ ξ F ′ = 0 (16)

k thnf / k f Pr ρ c p thnf / ρ c p f I ″ + f 0 I ′ + F θ 0 ′ ‐ 2 f 0 ′ I + θ 0 F ′ + β ξ I ρ c p thnf / ρ c p f = 0 , (17)

F ′ 0 = F 0 = 0 , I 0 = δ T I ′ 0 F ′ η → 0 , I η → 0 as η → ∞ (18)

In the next step, we follow the procedure of Harris et al. [40] and modify BCs by removing a condition from the original system. This action fixes the admissible ξ. The linearized stability equations are homogeneous and therefore always admit the trivial solution. For numerical analysis, a meaningful non-trivial solution is obtained by replacing one of BCs with a normalization condition that sets the perturbation magnitude at surface. Such normalization removes the trivial solution and allows the eigenvalue problem to be solved numerically by the shooting technique with the asymptotic BC used to find the associated ξ. It should be emphasized that the eigenvalue spectrum is independent of the normalization procedure, which only rescales the eigenfunction and does not affect the flow stability properties. Therefore, the condition F ″ 0 = 1 is used instead of the far-field condition F ′ η → 0 as η → ∞ to determine the smallest eigenvalue (ξ 1). The introduction of F ″ 0 = 1 is important in the stability analysis because it ensures that the perturbation solution is bounded and that unphysical exponential growth is avoided for large ξ. Thus, BC F ′ ∞ → 0 is replaced by the normalization condition F ″ 0 = 1. Any scalar multiple of F 0 η is a solution of the linearized equations as well, so we may set F ″ 0 = 1 without loss of generality. Therefore, Equations 16–18 can be considered and solved as an initial value problem.

4 Results and discussion

The validity of the current numerical method was verified through a comparison study using MATLAB’s bvp4c solver. Tables 2, 3 lists the values of (R e x 1 / 2 C f) obtained in the present investigation alongside those reported by Devi and Devi [40], Waini et al. [42] and Saleem et al. [41]. The excellent agreement between the datasets demonstrates the accuracy of the proposed Scheme. 1 is evident that the present results are in excellent agreement with the published data, and only negligible differences are due to numerical approximations. The close agreement of the present computations with the benchmark results confirms the correctness, stability, and high accuracy of the bvp4c implementation. Hence, the numerical method developed is reliable and can be used with confidence for further analysis of the physical problem.

TABLE 2. Comparison of computed R e C f values from the present approach with published results of [40, 42] for S = 0, λ = 1, and ϕ 3 = 0.

TABLE 3. Comparison of computed R e C f values from the present approach with published results of [41] for S = 0, λ = 1, and ϕ 1 = ϕ 3 = 0.

Figures 3, 4 depict the effects of λ, and the C u O particle loading ϕ 3, on the dimensionless R e x 1 / 2 C f, and the normalized heat-transfer rate R e x − 1 / 2 N u x, respectively. The results indicate that for a particular range of λ, two distinct solution branches emerge, when λ c i ≤ λ where i = 1 , 2 , 3, which implies that under the same physical conditions, two possible flow states can coexist, no solution exist when λ < λ c i where i = 1 , 2 , 3. Figure 3 indicates that λ has an important effect on the variation of R e x 1 / 2 C f. The C f x shows a significant change with the increase in λ (toward shrinking) and the increase in the volume fraction of ϕ 3 increases the magnitude of R e x 1 / 2 C f for both solution. The critical value of λ is also obtained beyond which the dual solutions tend to merge indicating a possible transition or loss of steady similarity solutions. The increased λ c i indicates that THNFs can maintain the BL flow throughout a large operating range, which is beneficial for industrial cooling and coating applications. Figure 4 also shows that R e x − 1 / 2 N u x is strongly affected by λ and ϕ 3. The heat transfer rate for first sol. Increases with increasing λ, while second sol. Has comparatively lower heat transfer characteristics. The existence of trihybrid nanoparticles increases the rate of heat transfer and the suggested fluid can be utilized in thermal management, electronic cooling and heat exchanger applications. Moreover, the increasing trend of R e x − 1 / 2 N u x with the surge of volume fraction of C u O nanoparticles ϕ 3 indicates the enhancement of thermal conductivity and heat transfer performance because of nanoparticles existence. The outcomes display that the λ effects and the nanoparticle concentration have a significant control over the overall flow resistance and thermal transport characteristics of the system.

Figures 5, 6 show the impact of S on the velocity (f ′ η) and temperature profiles (θ η) respectively. In Figure 5 we can see clearly that the velocity profile f ′ η decreases with increasing suction parameter for first sol. Physically, wall suction removes the low-momentum fluid accumulated in the near surface region, and thus suppresses momentum diffusion and decreases the momentum BL thickness. This increases the wall shear, first solution tends to reach the free-stream velocity faster and second solution exhibits the opposite behaviour because of its different flow structure related to the dual-solution regime. Moreover, the increased BL control by wall suction is useful for coating, drying, cooling and polymer extrusion processes. From Figure 6, findings reveal that θ η reduces for both solutions as the suction parameter increases. This means that the cooling is enhanced by extracting more thermal energy from BL. From a practical point of view, the reduction in thermal BL thickness due to wall suction improves the cooling performance; which is beneficial for heat exchangers, electronic cooling and thermal management systems. Figures 7, 8 depict the influence of the composition and concentration of nanoparticles on the f ′ η and θ η profiles. The results show that the trihybrid nanofluid (ϕ 1 = ϕ 2 = ϕ 3 = 0.1) comprising alumina, copper, and copper oxide nanoparticles exhibits the maximum velocity gradient near the wall and the minimum temperature profile among the cases considered. Physically, the synergistic interaction among the three nanoparticles increases both the effective viscosity and thermal conductivity, significantly affecting the thermophysical properties of the suspension. The higher viscosity alters the momentum diffusion, resulting in an increased velocity gradient near the wall. The enhanced thermal conductivity leads to greater thermal diffusion, producing a thinner thermal BL and increased heat transport from the surface. Thus, THNF displays superior momentum and heat transfer characteristics compared to HNF (alumina-copper, ϕ 1 = ϕ 2 = 0.1 , ϕ 3 = 0) and the mono nanofluid (alumina only, ϕ 1 = 0.1 , ϕ 2 = ϕ 3 = 0).

The effect of heat sink/source parameter (β) and thermal slip condition (δ T) on temperature profile θ η is shown in Figure 9. As β increases (i.e., heat source with β > 0 the θ η profile increases for both first and second sols., i.e., heat addition increases the fluid temperature in BL. Also, it is shown that the inclusion of a thermal slip condition (δ T = 0.1 leads to lower temperature profiles than the no-slip case (δ T = 0), indicating that thermal slip at the boundary reduces heat transfer to the fluid. The effect of thermal slip becomes more significant as the value of β increases. An examination of the results indicates that β and δ T are of significance in the determination of BL temperature distribution. From a practical point of view, the combined control of β and δ T provides an effective mechanism for controlling the surface temperature, which is useful for electronic cooling, thermal insulation and heat exchanger applications.

The linearized eigenvalue problem, Equations 16 and 17 subject to BC Equation 18, was solved using the MATLAB bvp4c solver with a convergence tolerance of 10 − 8, a computational domain of η ∞ = 20, and a starting mesh of 200 points. Figure 10 shows the smallest eigenvalues (ξ 1) for the two solutions as functions of the parameter λ, with all other parameters being constant. These results indicate that ξ 1 are positive in first region of sol (blue), so these solutions are stable. When we say stable in this context, we mean that any small perturbation to the system will decay in time, and the system will return to its equilibrium state. The second region of solution (second sol. red) on the contrary is characterized by negative eigenvalues, which means that it is unstable. Any small perturbation here will increase with time, driving the system away from equilibrium. The critical value λ c = − 1.3201 is the bifurcation value between stable and unstable branches. Henceforth, stability analysis confirms that only first sol. Is physically meaningful for steady-state operation, since second sol. Is intrinsically unstable due to the sign of its smallest eigenvalues.

4.1 GRA

In the previous section, the individual parameters have been investigated using numerical simulations, with respect to their effects on the velocity, temperature, skin friction coefficient and Nusselt number. In a real thermal system, however, several operating variables and several performance responses must be considered simultaneously when designing the system. Choosing the optimum combination of skin friction coefficient and heat transfer rate from the many numerical simulations is not easy since these two parameters are sometimes mutually conflicting. Thus, the numerical outputs of the stable solution branch are taken as the input data set for the Taguchi design, Grey Relational Analysis (GRA) and Response surface methodology (RSM). These methods allow optimization with multiple responses, allow calculation of the relative importance of the various governing parameters, and allow the establishment of predictive equations that can be used in engineering design without having to numerically re-solve the governing equations over and over again.

GRA, developed from Grey System Theory by Deng, is a powerful multi-response optimization technique employed when several quality characteristics must be optimized simultaneously. In manufacturing and engineering applications, different performance measures often exhibit conflicting behavior. Therefore, GRA converts multiple responses into a single performance index known as the Grey Relational Grade (GRG), enabling straightforward determination of optimal process parameters. Table 4 presents the parameters/factors involved in our study and their levels. We have taken five factors and four levels. In the light of these factors and levels our L16 orthogonal array of factors and the responses are shown in Table 5.

TABLE 4. Factors and their levels.

TABLE 5. L16 orthogonal array layout and performance results.

The procedure and algorithm of GRA for this study can be depicted in Figure 11. At first stage, we start from design of experiments (Tables 4, 5). Then we calculate the signal to noise ratio (S/N ratio) using the following formula:

The larger, the better S/N ratio is calculated from the following formula and values are presented in Table 5.

S / N ratio for N u x 0 = ‐ 10 × log 10 1 n * ∑ i = 1 n * 1 y ij 2 (19)

S / N ratio for C f x 0 = ‐ 10 × log 10 1 n * ∑ i = 1 n * y ij 2 (20)

Once, S/N ratios for N u x and C f x computed we normalized our data from Min/Max formula as follows:

x j * ω = y j ω − min ⁡ y j ω max ⁡ y j ω − min ⁡ y j ω (21)

x j * ω = max ⁡ y j ω − y j ω max ⁡ y j ω − min ⁡ y j ω (22)

Where ω = 2 are the number of response variables (C f x 0 , N u x 0 and j = 1 , 2 , 3 , … , 16 are the runs of the decision matrix.

After this step, we computed the grey relational coefficients with the aid of deviation sequences. The computational values are presented in Table 6 and formula is mentioned as below:

τ x j * ω , x ° * ω = Δ min ω + ξ Δ max ω Δ 0 j ω + Δ max ω (23)

TABLE 6. S/N ratio Normalization, deviation sequences, grey relational coefficients and grades.

Here, Δ 0 j ω = x ° * ω − x ω is the deviation sequence, the distinguishing coefficient ξ is taken 0.5.

Grey relational grade is calculated from Equation 24 and values are present in Table 6.

α j x ° , x j * = 1 n ∑ ω = 1 n τ x j * ω , x ° * ω (24)

Now, we rank the GRG grade. From Table 6, it is worthy to note that Run 11 is the best and ranked as 1 and Run 10 is the worst and ranked at 16^th place. Experiment 11 with Grey Relational Grade (GRG = 0.9345) is the optimal stable solution from the Grey Relational Analysis. The governing parameters having the optimal values are S = 2.6, λ = −0.5, δ t ​ = 0, β = 0, and ϕ = 0.30. This parameter combination is the best compromise which maximises the overall heat transfer performance of the trihybrid nanofluid system for stable solution branch with respect to the skin friction coefficient as well as Nusselt number. From this Table it is came to know that the parameter S has it optimum value at Level 4, for λ , δ t , β the optimum values are obtained from Level 1 and ϕ has its optimum value at Level 4. So, if we set these levels to parameters then we can achieve the optimum values of the response variables C f x 0 , N u x 0. From the results of GRA, it can be seen that there was a significant difference in the overall performance of the sixteen experimental runs. The best Grey Relational Grade (GRG = 0.9345) was obtained in the 11th experiment, which was thus ranked first and considered as the optimum operation condition from all the different parameter combinations tested. The 16th experiment produced the second highest GRG value of 0.7661, followed by the 6th experiment (GRG = 0.7159) and the 14th experiment (GRG = 0.7077) which were ranked third and fourth, respectively. The top ranked experiments possess very good multi-response features due to their normalized responses being closer to the ideal sequence, which leads to higher Grey Relational Coefficients, and consequently larger Grey Relational Grades. The lowest GRG value (0.3506) was obtained in the 10th experiment, which ranked 16th position with respect to the simultaneous optimization objectives. The minimum Grey Relational Grade (GRG = 0.3506) was obtained for Experiment 10, thus it was considered as the poorest performing one. The smaller the GRG the further the normalized responses are from the ideal reference sequence in Grey Relational Analysis. In this case, Experiment 10 was the least satisfactory compromise between the skin friction coefficient and Nusselt number, leading to the lowest multi-response optimization performance for the sixteen experimental runs. Likewise, the relatively poor overall performance of the 8th and 7th experiments, with corresponding GRG values of 0.3870 and 0.4166, respectively, demonstrated that these experiments had poorer performance. The ranking analysis shows that the 11th experimental condition offers the best compromise of the responses considered, while the lower ranked experiments are quite far from the ideal performance characteristics. Hence it can be concluded that the optimum combination of the parameters for the best overall performance of the system is the one corresponding to Experiment 11.

The mean of the grey relational grade for each level of parameters calculated in Table 7. The larger GRG values are bold in Table 7. Based on the parametric optimization and variance analysis, the most influential parameter in controlling the performance of the system is (ϕ = ϕ 3, which has a gradient change from maximum to minimum of 0.29588 with the first rank. The thermal BL thickness parameter δ t, is second with significant impact variation of 0.27339 at the tested levels. The parameters such as velocity slider coefficient λ and suction/injection parameter S show moderate effect on the Grey Relational Grade (GRG) and rank third and fourth respectively. The factor β shows the least significance with minimum variation response of 0.02795 and it is ranked at the last fifth position.

TABLE 7. Response table for grey relational grade.

4.2 Response surface analysis (RSM)

Table 8 presents the factors and their level. Low level of factors is denoted by − 1, moderate level of factors by 0 and high level of factors by 1. Corresponding to these levels, coded values of these factors and relative values of C f x and N u x are presented in Table 9. Table 10 presents the full ANOVA which is used to determine the statistically significance of the factors involved in response function C f x 0 . From Table 10 it is worthy to notice that the over all p value is 0.002 which indicates that the regression model is highly significant and adequate to describe the response function in term of factors. All factors who has p-value greater than 0.05 are considered to be insignificant and after imposing the said condition the regression equation of C f x 0 can be written as Equation 25. Individual and combined contribution of factors on N u x is depicted in Table 11. Overall p-value is less than 0.05 therefore model of N u x 0 is statistically significant. The response function N u x 0 can be regress from Equation 26.

C f x 0 = ‐ 135 + 103.2 λ + 103.6 ϕ + 25.1 S · λ ‐ 27.6 S · ϕ + 43.8 λ · ϕ + 40.9 δ t · β ‐ 38.9 β · ϕ (25)

TABLE 8. Factors and their levels.

TABLE 9. Coded values and responses.

TABLE 10. Analysis of variance for C f x

Regression equation in the light of above Table 10 can be written as.

TABLE 11. Analysis of variance for N u x.

The regression model showed very good predictive capacity with the values of R ² = 0.9149, adjusted R ² = 0.8603 and the predicted R ² = 0.8440, indicating that the model was adequate to model the skin friction coefficient, Cf_x (0). From Figure 12, detailed residual diagnostics were conducted to validate the model. The normal probability plot indicated that the residuals were approximately normally distributed with some minor deviations at the tails. The resulting histogram (which is the residual plot) was approximately symmetric about zero, meaning that the model does not have a systematic tendency to overpredict or underpredict the numeric answers from the MATLAB bvp4c simulations. In addition, the residuals against the fitted values showed no obvious systematic patterns or funnel-shaped graph indicating the validity of homoscedasticity and validity of the polynomial regression model. The residuals also were randomly scattered about zero as in the residuals versus observed order graph which indicates no sequential, computational or time-dependent bias. The results validate the developed regression model of prediction and optimization.

After computation, we get R 2 = 99.56 % , R 2 − a d j = 98.76 % , R 2 − p r e d = 92.05 % which shows that our model is adequate to describe the needs of the study. The residual diagnostic plots of this regression model for the local Nusselt number N u x 0, were validated using the normal probability plot, the histogram of the residuals, the residuals versus the fitted number and the residuals versus the order of observations (Figure 13). The normal probability plot of the residuals was consistent with the reference line, except for some minor deviations at the tails, meaning that the normality assumption was reasonably met. The prediction errors were small and unbiased, as shown by the histogram, which had a bell-shaped distribution around the zero with an approximately equal distribution on both sides of zero. Furthermore, there was no trend in the residuals plots against the fitted values, nor was there a funnel-shaped pattern, which again indicated constant error variance (homoscedasticity) and the second-order regression model was appropriate. The residuals against the order of the observations also showed random fluctuation around the zero value without any cyclic or sequential patterns, thus showing randomness and lack of computational and systematic bias. These diagnostics elucidate the validity and the predictive power of the developed regression model.

Regression equation of N u x 0 in the light of above Table 10 is:

N u x 0 = 13.1 ‐ 4.7 S ‐ 3.6 λ ‐ 13.4 δ t + 73.76 δ t 2 ‐ 19.35 S · δ t ‐ 9.67 λ · δ t + 14.86 δ t · β + 34.26 δ t · ϕ (26)

From the response surface analysis, the suction/injection parameter S and the stretching/shrinking parameter λ are seen to have a significant effect on the thermal and hydrodynamic parameters of the nanofluid flow (Figure 14a). There is a monotonic decrease in N u x for enhancement of the surface heat transfer with both the increase in either of the parameters (Figure 14b). For higher values of S and λ, where C f x is increasingly negative, indicating a higher shear stress at the walls. The almost planar response surfaces and contour plots indicate that the effects of the parameters are mainly linear and that the interactions between the parameters, in the investigated operating range, are relatively small. The results are useful to optimize permeable surface systems where there is a need to achieve a balance between heat transfer enhancement and frictional losses. The combined effect of the suction parameter S and the boundary parameter δ t on C f x 0 is shown in the three-dimensional response surface (Figure 15a). The systematic negative orientation suggests that the surface shear stress is a resistive force in the opposite direction to the primary fluid motion. This is a typical example of BL flows over a shrinking sheet flow configuration λ < 0 in which the backward flow momentum of the wall results in large velocity gradients and local deceleration flow near the boundary. Parallel evenly spaced diagonal contour projections are seen in the bottom plane. This topological configuration is useful in verifying that the flow resistance is mainly linearly additive over these parameter ranges for the individual factors, namely, S and δ t. From this profile, it is noticed that there are no strong quadratic curvature or non-linear coupling effects in this contour configuration within these parameter ranges. The critical maximum drag configuration is for sure isolated at the critical point of minimum slip and maximum suction δ t = 0.0 , S = 2.6. The three-dimensional response surface shows the simultaneous effects of the suction parameter S and the boundary parameter δ t on N u x as shown in Figure 15b. The numerical values of N u x 0 span from approximately − 1.5 to 4.5. Importantly, the surface enters the negative region over a very large parameter space (the dark blue area Figure 15b). When N u x is negative, heat transfer is reversed and the fluid absorbs heat from the solid surface instead of the heat dissipating into the fluid. This phenomenon is usually observed when internal heat production (heat source behavior, β > 0) or strong flow retardation occurs and the flow temperature near the wall becomes higher than the wall temperature. In contrast to the linear relationships seen in C f x characteristics, the value of N u x 0 is strongly non-linear and parabolic with the parameter δ t. As δ t increases from 0 to 0.2 the heat transfer rate decreases but does so initially to some minimum value before undergoing a significant rise. The effect of the suction parameter S and the heat source/sink parameter β on C f x 0 is shown in Figure 16a. The response surface changes continuously and decreases in color from bright yellow to deep blue as the suction parameter S is increased from 2.3 to 2.6. This is a gradual algebraic decrease of C f x, and thus an increase in the absolute value of the wall shear stress. The parameter β is allowed to vary from 0 (no heat generation) to 1.0 (strong heat source configuration) and the value of C f x does not change noticeably. The two effects of the suction parameter S and the heat source/sink parameter β on N u x are shown in Figure 16b. If the suction parameter is fixed, and the value of β is increased from 0.0 to 1.0, then N u x increases significantly, and the response surface moves from the blue region to the bright yellow band. From a thermodynamic point of view, with regard to this particular boundary configuration, an increase in the heat parameter changes the thermal profiles of the ambient fluid.

The simultaneous effect of ϕ = ϕ 3, and the suction parameter S, on C f x is investigated through the response surface (Figure 17a). When the suction parameter is kept constant, and ϕ = 0.0 = ϕ 3 to ϕ = 0.25 = ϕ 3, the absolute value of the wall C f x is seen to increase in an algebraic manner. Surface plot (Figure 17b) is shown to illustrate the simultaneous effects of N u x due to nanoparticle solid volume fraction ϕ and the wall suction parameter S. As in the earlier thermal analyses, the surface enters the negative region where it is located in the high suction and low nanoparticle concentration area (dark blue quadrant). Figure 18a represents the simultaneous effects of λ and δ t on C f x. A step increase in the parameter λ leads to a significant algebraic rise in the value of C f x 0, from deep blue to bright yellow. The 3-D response surface shows the simultaneous effect of the stretching/shrinking parameter and the boundary parameter on N u x is shown in Figure 18b. The combined effect of the parameter for heat source/sink β and that of the stretching/shrinking λ on C f x is shown in Figure 19a. All C f x values fall within the negative domain, indicating that the wall shear stress is always against the primary flow stream, as there is reverse flow at the shrinking wall λ < 0. A change in the parameter from −0.8 to −0.5results in an appreciable increase in C f x 0, towards less negative values since a lower shrinking rate leads to a less strong backward momentum transfer to the fluid and a reduction in the near-wall velocity gradient. The surface and the base contour projection of N u x are presented in Figure 19b to analyze the joint effects of the heat source/sink parameter β and the stretching/shrinking parameter λ. Figure 20a is used to investigate the simultaneous effects of ϕ = ϕ 3 and the stretching/shrinking parameter on C f x 0. The values for C f x are all within the negative range, which further confirms that the shear stress in the skin, is a continuous and opposing force to the main flow direction, as a result of the adverse momentum created by the shrinking boundary condition λ < 0. The three dimensional response surface plot along with the base projection plot of N u x with simultaneous effect of ϕ = ϕ 3 and stretching/shrinking parameter are plotted in Figure 20b. The N u x values are always positive over the entire domain, which indicates that convective heat transfer always removes energy from the sheet towards the surrounding nanofluid matrix. As ϕ = ϕ 3 is increasing from 0 to 0.2, the value of N u x 0 is found to be a consistent and monotonic increase as the addition of solid nanoparticles with a very high thermal conductivity increases the overall thermal conductivity of the base fluid.

Figure 21a presents C f x with respect to the boundary parameter δ t and the heat source/sink parameter β. The C f x values are all within the negative range of C f x, which indicates that the wall shear stress is in the opposite direction of the primary flow under the shrinking sheet mechanism. The N u x is studied as a function of the boundary parameter δ t and heat source/sink parameter β by the presented three dimensional response surface (Figure 21b). The C f x with respect to ϕ = ϕ 3 and the boundary parameter δ t are evaluated in Figure 22a. The range of C f x values is fully in the negative region which indicates in this case that the wall shear stress is a constant resisting force acting in the opposite direction of the primary flow as a result of the backward flow behavior of the shrinking wall condition. The simultaneous effects of ϕ = ϕ 3 and the boundary parameter on N u x 0 can be depicted from Figure 22b. The 3-D response surface is used to examine the interaction between the ϕ = ϕ 3 and β parameter on C f x is shown in Figure 23a. The C f x is found in the negative range, which is the domain of continuous opposing shear stress at the solid boundary, which is controlled by the shrinking sheet mechanism. The effect of ϕ = ϕ 3 and the heat source/sink parameter on N u x is depicted in the form of a 3D view of the response surface and its base projection plot in Figure 23b.

5 Conclusion

The numerical scheme validated by MATLAB’s bvp4c is in good agreement with the literature confirming the accuracy and stability. For some intervals of the stretching/shrinking parameter λ there are two solutions and below a critical value no solutions. The growth of Re x 1 / 2 C f and Re x − 1 / 2 Nu x with an increase in λ and CuO nanoparticle fraction (ϕ 3) enhances flow resistance and heat transfer. The suction (S) reduces the velocity in first sol. and increases it in second sol. while it decreases the temperature in both solutions. Under higher velocity gradients and lower temperatures than hybrid and single-phase fluids, THNF (ϕ 1 = ϕ 2 = ϕ 3 = 0.1) shows improved transport. The source of heat (β) increases the temperature, and the thermal slip (δ T) decreases it for larger β. The best-optimized set of the parameters is: S (Level 4), λ (Level 1), δ_t (Level 1), β (Level 1), ϕ = ϕ 3 (Level 4), obtained using GRA. This is the best state of compromise between hydrodynamic and thermal performance of the nanofluid flow system. The main factor that affects the overall system performance is ϕ = ϕ 3, and the thermal BL thickness parameter (δ_t) is the second most influential one. It can be seen that the heat source/sink parameter (β) has the lowest influence on the Grey Relational Grade, and thus the smallest contribution to multi-response optimization. The heat transfer performance is enhanced by the increase in nanoparticle concentration and the suction and stretching/shrinking parameters have a significant effect on the thermal and hydrodynamic characteristics of THNF flow. Overall, the present findings provide useful guidelines for the design and optimization of advanced thermal management systems, heat exchangers, electronic cooling and other industrial heat transfer applications using THNFs. The present mathematical model can be extended in the future to include non-Newtonian base fluids, variable thermophysical properties, magnetic field effects, entropy generation, and experimental validation, to increase its practical applicability.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Balasubramaniam S, Digital University Kerala, India

Reviewed by: Nasir Shehzad, Independent Researcher, Rawalpindi, Pakistan

Atul Ray, Madhav Institute of Technology & Science Gwalior, India