Spherical Hybrid Nanoparticles for Homann Stagnation-Point Flow in Porous Media via Homotopy Analysis Method

Abstract

Non-axisymmetric stagnant-point flows for flat plates in porous media containing spherical Cu-Al₂O₃-H₂O nanoparticles are studied using the homotopy analysis method (HAM). The governing equations are transformed into three coupled non-linear ordinary differential equations through similarity transformations. A large degree of freedom is provided by HAM when selecting auxiliary linear operators. By transforming nonlinear coupled ordinary differential equations with variable coefficients into linear ordinary differential equations with constant coefficients, nonlinear coupled ordinary differential equations can be solved. Over the entire domain, these equations can be solved approximately analytically. The analysis involves a discussion of the impact of many physical parameters generated in the proposed model. The results have shown that skin friction coefficients of Cf_x and Cf_y increase with volume fraction of hybrid nanofluid and the coefficient of permeability increasing. For the axisymmetric case of γ = 0, when volume fraction, φ, φ₁, φ₂ = 0, 5%, 10%, 20%, Cf_x = Cf_y = 1.33634, 1.51918, 1.73905, 2.33449, it can be found that the wall shear stress values increase by 13.68%, 30.14%, and 74.69%, respectively. In response to an increase in hybrid nanofluid volume fractions, local Nusselt numbers Nu_x increase. Nu_x decrease and change clearly with the coefficient of permeability increasing in the range of γ < 0; the values of Nu_x are less affected in the range of γ > 0.

1. Introduction

In the analysis of fluid dynamics, properties of heat transfer have always been the focus of research, especially in the field of engineering, where enhancing heat transfer is very important. Hybrid nanofluids have been proposed and verified using thermal productivity calculations based on numerical and experimental results. Mixed convective Al₂O₃-water nanofluids were studied experimentally on inclined copper tubes by Momin et al. [1]. Studies and discussions of how nanoparticle concentration and power supply affect laminar flow thermal fields have been conducted. In a two-stage procedure, Sahoo et al. [2] prepared (Al₂O₃-SiC-TiO₂) ternary hybrid nanofluids based on water. As a result of this analysis, Al₂O₃-SiC-TiO₂ nanoparticles were equally distributed across each volume fraction sample of ternary hybrid nanofluids. Using numerical simulations and experiments, Zufar et al. [3] studied heat pipes and found that those including hybrid nanofluids perform well thermally. Al₂O₃-CuO and SiO₂-CuO hybrid nanofluids were characterized with respect to their thermodynamic conductivity and viscosity. Hybrid nanofluid flows with coupling stresses on the tensile surface were examined by Saeed et al. [4] using the homotopy analysis method. Nasir et al. [5] considered the effects of nonlinear thermal radiation on nanofluids based on water, hybrid, and ternary hybrids on stretched sheets using the homotopy analysis method. Hybrid nanofluids were used experimentally as heat transfer mediums in collectors and storage systems for solar energy applications by Yasmin et al. [6].

It is possible to improve the convective heat transfer characteristics of processes by using hybrid nanofluids in porous media, such as geothermal, oil flow filtration, all types of heat exchangers, and so on. A variety of hybrid nanofluids for a variety of applications showed increased efficiency and energy savings. Various aspects of hybrid nanofluids, including their synthesis, thermophysics, and heat transfer, have recently been reviewed by Sarkar et al. [7] and Babu et al. [8]. Although hybrid nanofluids might transfer heat better with proper hybridization, problems with preparation and stability still need to be overcome. Kasaeian et al. [9] conducted a comprehensive review of nanofluids and porous media for improving a thermal system’s heat transfer with varying structures and flow states. Das et al. [10] thoroughly reviewed experiments and numerical research on different nanofluids, as well as the most recent advancements in the study of thermal conductivity. Additionally, nanoparticles with a variety of types and sizes, solid volume fractions, different types of basic fluids, temperatures, and various mechanisms are discussed, all of which influence thermal conductivity [11,12]. A number of models and media were proposed by Tiwar et al. [13] for study of the thermal parameters of solar collectors in presence of nanofluids. Fluids used for heat transfer that consist of nanofluids and hybrid nanofluids can enhance the efficiency of solar collector absorption tubes. Sangapatanam et al. [14] presented a review of hybrid nanofluids’ preparation, stability, and characteristics of thermophysics. Moreover, how to control the pH and ultrasonic intensity of hybrid nanofluids has been discussed. Research on the effects of nanofluid type, nanoparticle concentration, and nanofluid depth on solar cell performance was reviewed by Modi et al. [15]. How nanofluids can be applied to solar cells and various enhancement techniques have also been discussed [16]. As part of their experimental investigation, Jana et al. [17] explored the impact of hybrid nanoparticles on fluid heat transfer. Mechanisms for improving thermal conductivity and the effects of nanoparticles on stability were also discussed. According to Suresh et al. [18], hybrid nanofluids exhibited an increase in thermal conductivity with increasing concentrations, and they found that the increase in viscosity was much higher than of the increase in thermal conductivity. Takabi et al. [19] numerically studied the effects of hybrid nanofluids as working fluids on thermal properties of a shell and its temperature field. Devi et al. [20] studied flows and heat transfer phenomena using hybrid nanofluids under magnetic fields on permeable tensile sheets. Hybrid nanofluids that improve the performance of suspended nanoparticles were discussed by Nabil et al. [21]. Bibi et al. [22] investigated flows and heat transfer of Cu, Al₂O₃, and TiO₂-H₂O nanofluids in vertical planes of anisotropic permeable saturated porous media at high Rayleigh numbers. Hayat et al. [23] considered nanoparticle volume fraction, magnetohydrodynamics, and heat generation–absorption effects on hybrid nanofluid flows through a stretched surface, and also compared the heat transfer properties of conventional and novel hybrid nanofluids. Saeed et al. [24] studied the thermal characteristics of Darcy–Forchheimer hydromagnetic hybrid nanofluids flowing in permeable stretch cylinders, including Brownian motion and the thermophoresis effect. In a porous medium, Wahid et al. [25] numerically simulated flows of Marangoni hybrid nanofluids, demonstrating the importance of porous media and hybrid nanofluids. Using magnetic fields and heat sources/sinks, Othman et al. [26] investigated the flows and heat transfer behaviors of carbon nanotubes on permeable exponentially shrinking surfaces. Khan et al. [27] considered the radiation-mixed convection flows of hybrid nanofluids through vertical cylinders with porous heat sources and sinks. Rostami et al. [28] studied the flows of mixture-based dusty hybrid nanofluids on stretched plates in porous media under action of magnetic fields, and concluded that heat transfer efficiency could be increased by improving thermal and electrical conductivity. An investigation into the free convection of hybrid nanosuspensions in geometrically triangular porous vessels with internal heat generation, inclined Lorentz force, and radiation was conducted by Ahmed et al. [29]. A significant part of the heat transfer efficiency and fluid flow is determined by the heated and wavy walls within the chamber. Porous infinite plates were used by Chu et al. [30] to convey the thermal effects of hybrid nanofluids on different nanoparticles. Copper (Cu), titanium oxide (TiO₂), and aluminum oxide (Al₂O₃) with water-based liquids were used in the characterization of hybrid nanoparticles.

There are many references on stagnation point flow, some of which are given below. Ariel [31] studied the two-dimensional stagnation-point flow problem of non-Newtonian liquids. Mahapatra et al. [32,33] considered the steady stagnation flows of viscoelastic liquids on the stretched boundary, and obtained similar solutions of the stagnation flows of magnetic fluids. Weidman [34,35] investigated a modification of Homann’s axisymmetric stagnation-point flows, and compared it to a rigid plate so that it was non-axisymmetric. By applying the homotopy analysis method, Nawaz et al. [36] studied radially stretched thin plates causing laminar boundary layer flows of nanofluids, considering thermophoresis and Brownian motion. Azam et al. [37] numerically studied the unsteady magnetic fluid stagnation-point flows of Carreau nanofluids on the expansion/contraction cylinder under the action of nonlinear thermal radiation, taking into account the condition that the mass flux of nanoparticles at the boundary of the Buongiorno model is zero. Ahmed et al. [38] considered rotational stagnation flows of Maxwell nanofluids in porous rotating disks. Modified Buongiorno nanofluid models were used to investigate the effects of variable thermal conductivities and heat sources on heat transfer characteristics. On porous stretched/contracted plates, Kho et al. [39] studied two-dimensional Hermann stagnation-point flows and heat transfer of hybrid nanofluids. Numerical analysis of hybrid nanofluid flows in porous media in a non-axisymmetric stagnation region was conducted by Waini et al. [40].

A wide range of problems involving strong nonlinearity can be solved with the homotopy analysis method (HAM) [41]. The HAM method cannot rely on small or large physical parameters to produce results, in contrast with perturbation methods. Higher-order approximate series solutions can also be selected via HAM using its expression form [42,43]. It is also possible to use HAM to simply control the convergence of series solutions. With regard to analytic and semi-analytic methods for nonlinear partial differential equations, HAM demonstrates the unusual convergence characteristics of series solutions [44,45]. The homotopy analysis method was used by Xu et al. [46] to investigate the growth of velocity boundary layers caused by pulse initiation wedges. Analytical approximations of momentum boundary layers were obtained that matched Keller box numerical results. With arbitrary stretching velocity on a flat plate, You et al. [47] investigated non-similar boundary layer flows of second-order fluids. Based on the results of HAM, local skin friction coefficients, local Nusselt numbers, and boundary layer thicknesses were considered in detail [48]. According to Sardanyés et al. [49], homotopy solutions were obtained for a nonlinear model of cancer that described the interaction of tumor cells with healthy and immune cells. Mathematica software package BVPH 2.0 with HAM was used by Farooq et al. [50] to solve the magnetic hydrodynamics Falkner–Skan flow of nanofluids in a semi-infinite domain. In the stagnation region, Mustafa et al. [51] observed Carson fluid flowing to the stretched sheet. A homotopy analysis solution was obtained by analyzing the entire space domain. Heat transfer characteristics with viscous dissipation were also analyzed. In a porous medium facing an exponentially stretched porous plate, Ali et al. [52] presented homotopy analytical solutions for boundary layers and heat transfer flow. Researchers studied magnetohydrodynamic secondary nanofluid flows through biaxially stretched surfaces in the work of Ramzan et al. [53]. In unsaturated porous media, Patel et al. [54] used the homotopy analysis method to give the nonlinear one-dimensional solution to Boussinesq’s equation. A magneticohydrodynamic biovector Walter-B nanofluid flow caused by stretched thin plates was considered by Hayat et al. [55]. A combination of the melting parameter, radiation parameter, magnetic parameter, Brownian motion parameter, Prandtl number, Peclet number, Dufour number, and Soret number was investigated. Using a nonlinear reaction–diffusion system, Al-Qudaha et al. [56] proposed an optimal homotopy analysis algorithm. In this linearization algorithm, homotopy series solutions are constructed by using Taylor series approximation of nonlinear equations. The results showed that linearization improved homotopy series solution accuracy and convergence. Bottona et al. [57] mathematically analyzed the vertical movement of two different spherical non-evaporative droplets along the centerline. Taking into account Blasius viscous flow, magnetohydrodynamic flow, boundary layer flows resulting from free convection, and Von Kármán swirling viscous flow, a comparison of homotopy renormalization (HTR) and homotopy analysis (HAM) was made by Yang et al. [58]. Similarly, it was found that HAM approximation was significantly more accurate than HTR approximation. Magnetohydrodynamic flows and heat transfer analysis of Casson fluids on an exponentially shrinking thin plate were carried out by Liu et al. [59]. To solve the forced Duffing equation, Prof. Liao [60] used the non-perturbation method based on homotopy analysis, namely “Directly Defining Inverse Mapping” (MDDiM). With this algorithm, all small denominators were completely avoided, and multiple limit cycles of the forced Duffing equation with high nonlinearity could be obtained successfully.

Using the homotopy analysis method (HAM), this study approximates analytical solutions for hybrid nanofluids flowing through flat sheets at a Homann stationary-point. With the introduction in Section 1, Section 2 presents the mathematical solution and derivative using HAM. A discussion and graphic illustration of results is included in Section 3. Lastly, a summary of the main findings is included in Section 4.

2. Mathematical Model and Method

The Cu-Al₂O₃-H₂O hybrid nanofluid model for Hohmann stagnation-point flows in porous media is shown in Figure 1, where z is the normal direction, and xy is the plane in Cartesian coordinate systems. The governing equations are (Weidman [34,35]; Waini et al. [40]):

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

(1)

$$\overline{u}\frac{\partial \overline{u}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{u}}{\partial \overline{y}}+\overline{w}\frac{\partial \overline{u}}{\partial \overline{z}}={\overline{u}}_e\frac{d{\overline{u}}_e}{d\overline{x}}+\frac{{\mu }_{hnf}}{{\rho }_{hnf}}\frac{{\partial }^2\overline{u}}{\partial {\overline{z}}^2}-\frac{{\mu }_{hnf}}{K{\rho }_{hnf}}\left(\overline{u}-{\overline{u}}_e\right),$$

(2)

$$\overline{u}\frac{\partial \overline{v}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{v}}{\partial \overline{y}}+\overline{w}\frac{\partial \overline{v}}{\partial \overline{z}}={\overline{v}}_e\frac{d{\overline{v}}_e}{d\overline{x}}+\frac{{\mu }_{hnf}}{{\rho }_{hnf}}\frac{{\partial }^2\overline{v}}{\partial {\overline{z}}^2}-\frac{{\mu }_{hnf}}{\overline{K}{\rho }_{hnf}}\left(\overline{v}-{\overline{v}}_e\right),$$

(3)

$$\overline{u}\frac{\partial \overline{T}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{T}}{\partial \overline{y}}+\overline{w}\frac{\partial \overline{T}}{\partial \overline{z}}=\frac{{\overline{k}}_{hnf}}{{\left(\rho c_p\right)}_{hnf}}\frac{{\partial }^2\overline{T}}{\partial {\overline{z}}^2},$$

(4)

Assuming

$$\begin{array}{l}\overline{u}=0,\overline{v}=0,\overline{w}=0,\overline{T}=\begin{array}{cc}{\overline{T}}_w\left(x\right)& \begin{array}{cc}at& \overline{z}=\end{array}\end{array}0\\ \overline{u}\to {\overline{u}}_e\left(x\right),\overline{v}\to {\overline{v}}_e\left(x\right),\overline{w}\to {\overline{w}}_e\left(x\right),\begin{array}{cc}\overline{T}\to {\overline{T}}_{\infty }& \begin{array}{cc}at& \overline{z}\to \end{array}\end{array}\infty ,\end{array}$$

(5)

In the formula, $\overline{u},\overline{v},\overline{w}$ are the velocity components, the external flow velocities are ${\overline{u}}_e\left(\overline{x},\overline{y}\right)=\overline{x}\left(\overline{a}+\overline{b}\right)$, ${\overline{v}}_e\left(\overline{x},\overline{y}\right)=\overline{y}\left(\overline{a}-\overline{b}\right)$, ${\overline{w}}_e\left(\overline{z}\right)=-2\overline{a}\overline{z}$, where $\overline{a},\overline{b}$ are shear-strain rates. The surface temperature is ${\overline{T}}_w={\overline{T}}_{\infty }+{\overline{T}}_0\overline{x}$ (${\overline{T}}_{\infty }$ ambient temperature, ${\overline{T}}_0$ characteristic temperature), $\overline{K}$ represents the permeability of porous media. $\overline{k}$ means the conductivity in temperature, $\rho$ represents fluid density, μ is dynamic viscosity, and $c_p$ means specific heat capacity coefficient.

Table 1. Nanoparticles and H₂O thermophysical properties [40].

Type	$\mathit{\rho } ({\mathbf{kg}/\mathbf{m}}^3)$	k (W/mK)	${\mathit{c}}_{\mathit{p}} (\mathbf{J}/\mathbf{kgK})$
H2O	997.1	0.613	4179
Cu	8933	401	385
Al2O3	3970	40	765

shows thermophysical properties of pure water, Cu, and Al₂O₃ nanoparticles.

Table 2. Applied model for thermophysical properties of spherical hybrid nanoparticles [40].

Property	Mathematical Relations
Density	${\rho }_{hnf}={\varphi }_1{\rho }_{n1}+{\varphi }_2{\rho }_{n2}+\left(1-\varphi \right){\rho }_f$
Heat capacity	${\left(\rho c_p\right)}_{hnf}={\varphi }_1{\left(\rho c_p\right)}_{n1}+{\varphi }_2{\left(\rho c_p\right)}_{n2}+\left(1-\varphi \right){\left(\rho c_p\right)}_f$
Dynamic viscosity	${\mu }_{hnf}=\frac{\mu }{{\left(1-\varphi \right)}^{2.5}}$
Thermal conductivity	$\frac{k_{hnf}}{k_f}=\frac{\frac{{\varphi }_1{k}_{n1}+{\varphi }_2{k}_{n2}}{\varphi }+2k_f+2\left({\varphi }_1{k}_{n1}+{\varphi }_2{k}_{n2}\right)-2\varphi k_f}{\frac{{\varphi }_1{k}_{n1}+{\varphi }_2{k}_{n2}}{\varphi }+2k_f-\left({\varphi }_1{k}_{n1}+{\varphi }_2{k}_{n2}\right)+\varphi k_f}$

provides the applied model for the thermophysical properties of spherical hybrid nanoparticles. The symbols ${\varphi }_1$ and ${\varphi }_2$ represent Cu and Al₂O₃ nanoparticles, respectively, and $\varphi ={\varphi }_1+{\varphi }_2$, where solid compositions are represented by subscripts $n1$ and $n2$.

Using the similarity transformation (Weidman [34,35]; Waini et al. [40]):

$$\begin{array}{l}\overline{u}=\overline{x}{f}^{\prime }\left(\eta \right)\left(\overline{a}+\overline{b}\right),\overline{v}=\overline{y}{g}^{\prime }\left(\eta \right)\left(\overline{a}-\overline{b}\right),\overline{w}=-\sqrt{\overline{a}{\nu }_f}\left. [f\left(\eta \right)\left(\overline{a}+\overline{b}\right)+g\left(\eta \right)\left(\overline{a}-\overline{b}\right)\right]\\ \theta =\frac{\overline{T}-{\overline{T}}_{\infty }}{{\overline{T}}_w-{\overline{T}}_{\infty }},\eta =\overline{z}\sqrt{\frac{\overline{a}}{{\nu }_f}}\end{array},$$

(6)

In the formula, ‘ means the derivative of η.

Substitute Equation (6) into Equations (1)–(4), and one has:

$$A_1{f}^{‴}+\left(1+\gamma \right)\left(f{f}^{″}+1-{f^{\prime }}^2\right)+\left(1-\gamma \right)g{f}^{″}-A_1\sigma \left(f^{\prime }-1\right)=0,$$

(7)

$$A_1{g}^{‴}+\left(1-\gamma \right)\left(g{g}^{″}+1-{g^{\prime }}^2\right)+\left(1+\gamma \right)f{g}^{″}-A_1\sigma \left(g^{\prime }-1\right)=0,$$

(8)

$$\frac{A_2}{\mathrm{Pr}}{\theta }^{″}+\left(1+\gamma \right)\left(f{\theta }^{\prime }-f^{\prime }\theta \right)+\left(1-\gamma \right)g{\theta }^{\prime }=0,$$

(9)

depending on

$$\begin{array}{l}f\left(0\right)=0,f^{\prime }\left(0\right)=0,g\left(0\right)=0,g^{\prime }\left(0\right)=0,\theta \left(0\right)=1\\ f^{\prime }\left(\eta \right)\to 1,g^{\prime }\left(\eta \right)\to 1,\theta \left(\eta \right)\to 0\begin{array}{c}\begin{array}{cc}at& \eta \to \end{array}\end{array}\infty ,\end{array}$$

(10)

In the formula, the reflective symmetries are obtained as $f\left(\eta ,\gamma \right)=g\left(\eta ,-\gamma \right)$ or $f\left(\eta ,-\gamma \right)=g\left(\eta ,\gamma \right)$, where γ is the ratio of strain shear rate. The coefficients of $A_1,A_2$, $\sigma$, and Pr are given by

$$A_1=\frac{{\mu }_{hnf}/{\mu }_f}{{\rho }_{hnf}/{\rho }_f},A_2=\frac{k_{hnf}/k_f}{{\left(\rho c_p\right)}_{hnf}/{\left(\rho c_p\right)}_f},\gamma =\frac{b}{a},\sigma =\frac{{\nu }_f}{aK},\mathrm{Pr}=\frac{{\mu }_f{\left(c_p\right)}_f}{k_f}.$$

(11)

The skin friction coefficients of $C{f}_x,C{f}_y$ and the local Nusselt numbers $N{u}_x$ are

$${\overline{C}}_{fx}=\frac{{\mu }_{hnf}}{{\rho }_f{u}_e^2}{\left(\frac{\partial \overline{u}}{\partial \overline{z}}\right)}_{\overline{z}=0},{\overline{C}}_{fy}=\frac{{\mu }_{hnf}}{{\rho }_f{v}_e^2}{\left(\frac{\partial \overline{v}}{\partial \overline{z}}\right)}_{\overline{z}=0},{\overline{Nu}}_x=-\frac{\overline{x}{k}_{hnf}}{k_f\left({\overline{T}}_w-{\overline{T}}_{\infty }\right)}{\left(\frac{\partial \overline{T}}{\partial \overline{z}}\right)}_{\overline{z}=0},$$

(12)

$$\begin{array}{l}{C}_{fx}=\sqrt{{\mathrm{Re}}_x\left(1+\gamma \right)}{\overline{C}}_{fx}=\frac{{\mu }_{hnf}}{{\mu }_f}{f}^{″}\left(0\right),C_{fy}=\sqrt{{\mathrm{Re}}_y\left(1-\gamma \right)}{\overline{C}}_{fy}=\frac{{\mu }_{hnf}}{{\mu }_f}{g}^{″}\left(0\right),\\ N{u}_x=\sqrt{\frac{1+\gamma }{{\mathrm{Re}}_x}}{\overline{Nu}}_x=-\frac{k_{hnf}}{k_f}{\theta }^{″}\left(0\right)\end{array},$$

(13)

In the formula, the local Reynolds numbers are ${\mathrm{Re}}_x=u_{e}x/{\nu }_f$, ${\mathrm{Re}}_y=v_{e}y/{\nu }_f$.

With HAM, it is possible to simplify nonlinear problems into linear ones in an infinite number of ways. It is generally accepted that boundary layer flows degenerate exponentially after infinity [61,62]. Based on Equations (13)–(14), $f\left(\eta \right),\theta \left(\eta \right)$ could be expressed by

$$f(\eta )={\displaystyle \sum _{t=0}^{+\infty }{\displaystyle \sum _{s=0}^{+\infty }{a}_{t,s}{\eta }^t{e}^{-\beta s\eta }},g(\eta )={\displaystyle \sum _{t=0}^{+\infty }{\displaystyle \sum _{s=0}^{+\infty }{b}_{t,s}{\eta }^t{e}^{-\beta s\eta }},\theta (\eta )={\displaystyle \sum _{t=0}^{+\infty }{\displaystyle \sum _{s=0}^{+\infty }{c}_{t,s}{\eta }^t{e}^{-\beta s\eta }}}}},$$

(14)

In the formula, $a_{t,s},b_{t,s},c_{t,s}$ are constant coefficients determined by HAM.

Equation (14) provides a convenient way to select initial guesses using Equation (10)

$$f_0(\eta )=\eta -\frac{1-e^{-\beta \eta }}{\beta },g_0(\eta )=\eta -\frac{1-e^{-\beta \eta }}{\beta },{\theta }_0(\eta )=e^{-\beta \eta },$$

(15)

with auxiliary linear operators

$$L_f=\frac{{\partial }^{3}F}{\partial {\eta }^3}-{\beta }^2\frac{\partial F}{\partial \eta },L_g=\frac{{\partial }^{3}G}{\partial {\eta }^{3G}}-{\beta }^2\frac{\partial G}{\partial \eta },L_{\theta }=\frac{{\partial }^2\Theta }{\partial {\eta }^2}-{\beta }^2\Theta ,$$

(16)

having the following properties:

$$\begin{array}{l}{L}_f\left. [C_0+C_1{e}^{-\beta \eta }+C_2{e}^{\beta \eta }\right]=0\\ L_g\left. [C_3+C_4{e}^{-\beta \eta }+C_5{e}^{\beta \eta }\right]=0,\\ L_{\theta }\left. [C_6\exp e^{-\beta \eta }+C_7{e}^{\beta \eta }\right]=0\end{array}$$

(17)

In the formula, C₀–C₇ are integral coefficients.

HAM deformation equations are constructed:

$$\begin{array}{l}\left(1-q\right)L_f\left. [F\left(\eta ;q\right)-f_0\left(\eta \right)\right]=q{h}_f{N}_f\left. [F\left(\eta ;q\right)\right]\\ \left(1-q\right)L_g\left. [G\left(\eta ;q\right)-g_0\left(\eta \right)\right]=q{h}_g{N}_g\left. [G\left(\eta ;q\right)\right],\\ \left(1-q\right)L_{\theta }\left. [\Theta \left(\eta ;q\right)-{\theta }_0\left(\eta \right)\right]=q{h}_{\theta }{N}_{\theta }\left. [\Theta \left(\eta ;q\right)\right]\end{array}$$

(18)

and

$$\begin{array}{l}F\left(0;q\right)=0,F^{\prime }\left(0;q\right)=0,F^{\prime }\left(\infty ;q\right)=1\\ G\left(0;q\right)=0,G^{\prime }\left(0;q\right)=0,G^{\prime }\left(\infty ;q\right)=1,\\ \Theta \left(0;q\right)=1,{\Theta }^{\prime }\left(\infty ;q\right)=0\end{array}$$

(19)

In the formula embedding parameter $q\in [0,1]$, the nonlinear operator based on governing Equations (7)–(9) is $N_f,N_g,N_{\theta }$. Using HAM, approximate analytically solutions can be fully determined:

$$\begin{array}{l}F\left(\eta ;0\right)=f_0\left(\eta \right),\hspace{1em}\begin{array}{c}F\left(\eta ;1\right)=f\left(\eta \right)\end{array}\\ G\left(\eta ;0\right)=g_0\left(\eta \right),\hspace{1em}\begin{array}{c}G\left(\eta ;1\right)=g\left(\eta \right).\end{array}\\ \Theta \left(\eta ;0\right)={\theta }_0\left(\eta \right),\hspace{1em}\begin{array}{c}\Theta \left(\eta ;1\right)=\theta \left(\eta \right)\end{array}\end{array}$$

(20)

Q contributes to mapping; as q changes from 0 to 1, mapping ensures that $F\left(\eta ;q\right),G\left(\eta ;q\right),\Theta \left(\eta ;q\right)$ deform continuously to exact solutions $f\left(\eta \right),g\left(\eta \right),\theta \left(\eta \right)$ from initial guesses $f_0\left(\eta \right),g_0\left(\eta \right),{\theta }_0\left(\eta \right)$. By Taylor’s theorem, $F\left(\eta ;q\right),G\left(\eta ;q\right),\Theta \left(\eta ;q\right)$ for power series expansion are:

$$\begin{array}{l}F\left(\eta ;q\right)=F\left(\eta ;0\right)+{\displaystyle \sum _{s=1}^{+\infty }{f}_s\left(\eta \right)}{q}^s\\ G\left(\eta ;q\right)=G\left(\eta ;0\right)+{\displaystyle \sum _{s=1}^{+\infty }{g}_s\left(\eta \right)}{q}^s,\\ \Theta \left(\eta ;q\right)=\Theta \left(\eta ;0\right)+{\displaystyle \sum _{s=1}^{+\infty }{\theta }_s\left(\eta \right)}{q}^s\end{array}$$

(21)

In the formula

$$f_s\left(\eta \right)={\left.\frac{1}{s!}\frac{{\partial }^{s}G\left(\eta ;q\right)}{\partial q^s}\right|}_{q=0},g_s\left(\eta \right)={\left.\frac{1}{s!}\frac{{\partial }^{s}G\left(\eta ;q\right)}{\partial q^s}\right|}_{q=0},{\theta }_s\left(\eta \right)={\left.\frac{1}{s!}\frac{{\partial }^{s}G\left(\eta ;q\right)}{\partial q^s}\right|}_{q=0}.$$

(22)

Series (21) can be substituted into zero order deformation Equation (18) and boundary Condition (19) as follows. In order to obtain m order deformation equation, equal power coefficients of q are used:

$$\begin{array}{l}{L}_f\left. [g\left(\eta \right)-{\chi }_m{f}_{m-1}\left(\eta \right)\right]=h_f{R}_{fm}\left(\eta \right), \begin{array}{c}m\ge 1\end{array}\\ L_g\left. [g\left(\eta \right)-{\chi }_m{g}_{m-1}\left(\eta \right)\right]=h_g{R}_{gm}\left(\eta \right), \begin{array}{c}m\ge 1,\end{array}\\ L_{\theta }\left. [g\left(\eta \right)-{\chi }_m{\theta }_{m-1}\left(\eta \right)\right]=h_{\theta }{R}_{\theta m}\left(\eta \right), \begin{array}{c}m\ge 1\end{array}\end{array}$$

(23)

subject to

$$\begin{array}{l}{f}_m\left(0\right)=0,\hspace{1em}\begin{array}{c}f{\prime }_m\left(0\right)=0, \begin{array}{c}f{\prime }_m\left(\eta \right)=0\end{array}\end{array}\\ g_m\left(0\right)=0,\hspace{1em}\begin{array}{c}g{\prime }_m\left(0\right)=0, \begin{array}{c}g{\prime }_m\left(\eta \right)=0,\end{array}\end{array}\\ {\theta }_m\left(0\right)=1,\hspace{1em}\begin{array}{c}\theta {\prime }_m\left(\infty \right)=0\end{array}\end{array}$$

(24)

In the formula

$$\begin{array}{ll}{R}_{fm}\left(\eta \right)=& A_1\frac{{\partial }^3{f}_{m-1}}{\partial {\eta }^3}+\left(1+\gamma \right)\left({\displaystyle \sum _{n=1}^{m-1}{f}_{m-1-n}}\frac{{\partial }^2{f}_{m-1}}{\partial {\eta }^2}-{\displaystyle \sum _{n=1}^{m-1}\frac{\partial f_{m-1-n}}{\partial \eta }}\frac{\partial f_{m-1}}{\partial \eta }\right)\\ & +\left(1-\gamma \right){\displaystyle \sum _{n=1}^{m-1}{g}_{m-1-n}}\frac{{\partial }^2{f}_{m-1}}{\partial {\eta }^2}-A_1\sigma \frac{\partial f_{m-1}}{\partial \eta }+\left(1-C_m\right)\left(1+\gamma +A_1\sigma \right),\end{array}$$

(25)

$$\begin{array}{ll}{R}_{gm}\left(\eta \right)=& A_1\frac{{\partial }^3{g}_{m-1}}{\partial {\eta }^3}+\left(1-\gamma \right)\left({\displaystyle \sum _{n=1}^{m-1}{g}_{m-1-n}}\frac{{\partial }^2{g}_{m-1}}{\partial {\eta }^2}-{\displaystyle \sum _{n=1}^{m-1}\frac{\partial g_{m-1-n}}{\partial \eta }}\frac{\partial g_{m-1}}{\partial \eta }\right)\\ & +\left(1+\gamma \right){\displaystyle \sum _{n=1}^{m-1}{f}_{m-1-n}}\frac{{\partial }^2{g}_{m-1}}{\partial {\eta }^2}-A_1\sigma \frac{\partial g_{m-1}}{\partial \eta }+\left(1-C_m\right)\left(1-\gamma +A_1\sigma \right),\end{array}$$

(26)

$$\begin{array}{ll}{R}_{\theta m}\left(\eta \right)=& \frac{A_2}{\mathrm{Pr}}\frac{{\partial }^2{\theta }_{m-1}}{\partial {\eta }^2}+\left(1+\gamma \right)\left({\displaystyle \sum _{n=1}^{m-1}{f}_{m-1-n}}\frac{\partial {\theta }_{m-1}}{\partial \eta }-{\displaystyle \sum _{n=1}^{m-1}{\theta }_{m-1}\frac{\partial f_{m-1-n}}{\partial \eta }}\right)\\ & +\left(1-\gamma \right){\displaystyle \sum _{n=1}^{m-1}\frac{\partial {\theta }_{m-1}}{\partial \eta }{g}_{m-1-n}},\end{array}$$

(27)

$$C_m=\left\{\begin{array}{l}\begin{array}{cc}1& m>1\end{array}\\ \begin{array}{cc}0& m=0.\end{array}\end{array}\right.$$

(28)

In Equations (25)–(27), the right-hand side is given by Equations (7)–(9), which are actually an ordinary differential equation about η. The particular solutions of Equations (25)–(27) are as follows:

$$f_m^{*}\left(\eta \right)={\chi }_m{f}_{m-1}\left(\eta \right)+\hslash {\displaystyle {\int }_0^{\eta }\left\{{\displaystyle {\int }_0^{\eta }\left. [e^{-\beta \eta }{R}_{fm}\left(\eta \right)e^{\beta s}ds\right]d\eta }\right\}}d\eta .$$

(29)

$$g_m^{*}\left(\eta \right)={\chi }_m{g}_{m-1}\left(\eta \right)+\hslash {\displaystyle {\int }_0^{\eta }\left\{{\displaystyle {\int }_0^{\eta }\left. [e^{-\beta \eta }{R}_{gm}\left(\eta \right)e^{\beta s}ds\right]d\eta }\right\}}d\eta .$$

(30)

$${\theta }_m^{*}\left(\eta \right)={\chi }_m{\theta }_{m-1}\left(\eta \right)+\hslash {\displaystyle {\int }_0^{\eta }\left. [e^{-\beta \eta }{R}_{\theta m}\left(\eta \right)e^{\beta s}ds\right]d\eta }.$$

(31)

Equation (17) leads to the following general solution:

$$f_m\left(\eta \right)=f_m^{*}\left(\eta \right)+C_{0,m}+C_{1,m}\exp \left(-\beta \eta \right)+C_{2,m}\exp \left(\beta \eta \right).$$

(32)

$$g_m\left(\eta \right)=g_m^{*}\left(\eta \right)+C_{3,m}+C_{4,m}\exp \left(-\beta \eta \right)+C_{5,m}\exp \left(\beta \eta \right).$$

(33)

$${\theta }_m\left(\eta \right)={\theta }_m^{*}\left(\eta \right)+C_{6,m}\exp \left(-\beta \eta \right)+C_{7,m}\exp \left(\beta \eta \right).$$

(34)

Consequently, an infinite number of linear ordinary differential Equation (23) have constant coefficients derived from the original nonlinear partial differential Equations (7)–(9). The solution of nonlinear partial differential equations involving variable coefficients is undoubtedly easier than solving linear ordinary differential equations involving constant coefficients. There is no requirement for any small or large physical parameters in this transformation, unlike the perturbation method. There is no obvious relationship between the linear terms in the original Equations (7)–(9) and the auxiliary linear Operator (16). The main reason for this is that the homotopy analytic method is very flexible in terms of auxiliary linear operators, unlike other analytical methods. The linear ordinary differential Equation (23) with variable coefficients cannot be converted into the nonlinear coupled ordinary differential Equations (7)–(9) without this degree of freedom. The following series solutions of Equations (7)–(9) can be obtained.

$$f\left(\eta \right)=f_0\left(\eta \right)+{\displaystyle \sum _{m=1}^{+\infty }{f}_m\left(\eta \right)},g\left(\eta \right)=g_0\left(\eta \right)+{\displaystyle \sum _{m=1}^{+\infty }{g}_m\left(\eta \right)},\theta \left(\eta \right)={\theta }_0\left(\eta \right)+{\displaystyle \sum _{m=1}^{+\infty }{\theta }_m\left(\eta \right)}.$$

(35)

noting that $f\left(\eta \right),g\left(\eta \right),\theta \left(\eta \right)$ given by the homotopy analysis method contain three unknown convergence control parameters $h_f,h_g,h_{\theta }$. The series solution relies on these to ensure convergence. The mean residual errors of the kth order approximations are defined as follows:

$${\epsilon }_k^f\left(h_f,h_g,h_{\theta }\right)=\frac{1}{N+1}{{\displaystyle \sum _{i=0}^N\left. [{\left.N_f\left({\displaystyle \sum _{m=0}^k{f}_m}\right)\right|}_{\eta =i\beta \eta }\right]}}^2,$$

(36)

$${\epsilon }_k^g\left(h_f,h_g,h_{\theta }\right)=\frac{1}{N+1}{{\displaystyle \sum _{i=0}^N\left. [{\left.N_g\left({\displaystyle \sum _{m=0}^k{f}_m},{\displaystyle \sum _{m=0}^k{g}_m}\right)\right|}_{\eta =i\beta \eta }\right]}}^2,$$

(37)

$${\epsilon }_k^{\theta }\left(h_f,h_g,h_{\theta }\right)=\frac{1}{N+1}{{\displaystyle \sum _{i=0}^N\left. [{\left.N_{\theta }\left({\displaystyle \sum _{m=0}^k{f}_m},{\displaystyle \sum _{m=0}^k{g}_m},{\displaystyle \sum _{m=0}^k{\theta }_m}\right)\right|}_{\eta =i\beta \eta }\right]}}^2,$$

(38)

For the original governing Equations (7)–(9), an approximation of the kth order has a total error defined as follows:

$${\epsilon }_k^{tol}\left(h_f,h_g,h_{\theta }\right)={\epsilon }_k^f\left(h_f,h_g,h_{\theta }\right)+{\epsilon }_k^g\left(h_f,h_g,h_{\theta }\right)+{\epsilon }_k^{\theta }\left(h_f,h_g,h_{\theta }\right).$$

(39)

As a result of the kth order approximation, the optimal value for $h_f,h_g,h_{\theta }$ can be determined as the minimum of total error of ${\epsilon }_k^{tol}$.

3. Results Analysis and Discussion

Based on the HAM results, Cu-Al₂O₃-H₂O spherical hybrid nanoparticles for Homann stagnation-point flows of flat sheets in porous media are solved approximately and analytically. The convergence control parameter h and auxiliary parameter β in the series Solution (35) ensure the convergence of order solutions. By examining the curves of the control parameters, we are able to determine the optimal range for the auxiliary control parameter, which is crucial to ensuring that the series solution converges. Therefore, if the series Solution (35) is converging, it follows that the derivative $f_{\eta \eta }\left(\eta \right),g_{\eta \eta }\left(\eta \right),{\theta }_{\eta }\left(\eta \right)$ must also be converging. In order to simplify things, consider the convergence of the $f_{\eta \eta }\left(\eta \right),g_{\eta \eta }\left(\eta \right),{\theta }_{\eta }\left(\eta \right)$ series solutions. Because the original nonlinear partial differential Equations (7)–(9) are coupled, β must be greater than 2. Unless β and η are fixed, series solutions of $f_{\eta \eta }\left(\eta \right),g_{\eta \eta }\left(\eta \right),{\theta }_{\eta }\left(\eta \right)$ are only dependent on control parameter h. As shown in Figure 2, let $\beta =5$, and $f_{\eta \eta }\left(\eta \right),g_{\eta \eta }\left(\eta \right),{\theta }_{\eta }\left(\eta \right)$ are just power series of h whose convergence is dependent upon h when η = 0. There are two points that need to be emphasized. A first benefit of the HAM method is that it is a superior method to other analytical and semi-analytic methods. Nonlinear problems with strong nonlinearity can be converged accurately with HAM. In spite of large perturbations, HAM maintains independence from small or large physical parameters and presents convenient control over homotopy series solutions. A second advantage is that it converges rapidly when the convergence control parameters are optimized. With an increasing order of approximation, the residual total decreases for each case.

The analysis and discussion of results obtained using HAM is as follows. The analysis involves a discussion of the impact of many physical parameters generated in the proposed model. When γ = 0, $\varphi ={\varphi }_1={\varphi }_2=0$, and σ = 0, it can be found that $f\left(\eta \right)=g\left(\eta \right)$ represents the axisymmetric Homann stagnation-point flow. When γ = 0 (axisymmetric), $\varphi ={\varphi }_1={\varphi }_2=0$ (regular fluid), and $\sigma =0$ (non-porous medium), $f^{″}\left(0\right)$ = 1.311608, compared with the $f^{″}\left(0\right)$ = 1.311938 of Waini et al. [40], Wang [63], and Soid et al. [64], and the relative error is not more than 0.025%.

Table 3. Skin friction coefficients $C{f}_x,C{f}_y$ and Nusselt numbers $N{u}_x$ with various values of γ, $\varphi$ and σ when Pr = 6.2.

γ	$\mathit{\phi }$	σ	$\mathit{C}{\mathit{f}}_{\mathit{x}}$ (Refs. [40,62,63])	HAM 20th	$\mathit{C}{\mathit{f}}_{\mathit{y}}$ (Refs. [40,62,63])	HAM 20th	$\mathit{N}{\mathit{u}}_{\mathit{x}}$ (Refs. [40,62,63])	HAM 20th
0	0	0	1.311938	1.311608	1.311938	1.311608	1.806069	1.810147
5			3.038940	3.036096	−0.894909	−0.902242	3.938146	3.998352
−5			−0.894909	−0.902242	3.038940	3.036096	3.074275	3.084240
−5	2%		−0.966699	−0.966456	3.282727	3.281901	3.203682	3.210915
	4%		−1.039700	−1.039438	3.530622	3.529834	3.330939	3.669205
	2%	2	−0.027231	−0.027224	3.569397	3.568499	2.788727	2.795023

presents a comparison of corresponding values of $f_{\eta \eta }\left(0\right)$ and ${\theta }_{\eta }\left(0\right)$ with γ, $\varphi$, σ. A 20th order homotopy solution is adopted without higher order calculation, and the convergence result is accurate enough, especially for the large values of γ, $\varphi$, and σ. As shown in Table 3, skin friction coefficients $C{f}_x,C{f}_y$ and Nusselt number $N{u}_x$ with various values of γ, $\varphi$, σ when Pr = 6.2 are calculated and compared with the results of Waini et al. [40]. The consistency of results is very good. When Pr = 6.2, σ = 0, and the volume fraction $\varphi ={\varphi }_1={\varphi }_2=0$, the velocities of $f^{\prime }\left(\eta \right),g^{\prime }\left(\eta \right)$ and the temperature $\theta \left(\eta \right)$ distributions under different values of shear-strain rate ratios γ = 0, ±1, ±3 are presented in Figure 2. The behaviors of flow fields by changing the shear–strain rate ratios γ are studied. When $\gamma >0$, $f^{\prime }\left(\eta \right),g^{\prime }\left(\eta \right)$ increase with an increase in γ; $\theta \left(\eta \right)$ decreases with an increase in γ. When $\gamma <0$, $f^{\prime }\left(\eta \right),g^{\prime }\left(\eta \right)$ decrease as γ decreases; $\theta \left(\eta \right)$ decreases with the decrease in γ generally. $f^{\prime }\left(\eta \right)$ flows in the reverse direction towards the plate at $\gamma =-3$, or $g^{\prime }\left(\eta \right)$ flows in the reverse direction towards the plate at γ = 3. Flows are inward near the stagnation zone, as shown in Figure 3. When Pr = 6.2 and $\varphi ={\varphi }_1={\varphi }_2=5\%$, the velocities of $f^{\prime }\left(\eta \right),g^{\prime }\left(\eta \right)$ and the temperature $\theta \left(\eta \right)$ curves with different γ and σ values are shown in Figure 4a,b. $f^{\prime }\left(\eta \right),g^{\prime }\left(\eta \right)$, $\theta \left(\eta \right)$ increase as σ increases. When the volume fraction $\varphi ,{\varphi }_1,{\varphi }_2$ of spherical hybrid nanoparticles Cu-Al₂O₃-H₂O increases, the coefficients of skin friction and local Nusselt numbers increase. The Prandtl number is fixed at 6.2, and does not consider the variations.

Figure 5a,b illustrates variations in the skin friction coefficients $C{f}_x,C{f}_y$ with γ in the range of $-6\le \gamma \le 6$ for various $\varphi ={\varphi }_1={\varphi }_2=0$, 5%, 10%, 20%, and σ = 0, 1, 2, 5 when Pr = 6.2. $C{f}_x,C{f}_y$ show a symmetric pattern, where the line of symmetry lies at γ = 0 in the axisymmetric case. When γ = 0, σ = 0, and $\varphi ={\varphi }_1={\varphi }_2=0$, $C{f}_x=f^{″}\left(0\right)$ = 1.311608 and $C{f}_y=g^{″}\left(0\right)$ = 1.311608, also as shown in Table 3. As shown in Figure 2a, $C{f}_x=C{f}_y$ = 0, and the shear stress value is 0 when $\gamma =\pm 2.50319$. When $\gamma <-2.503193$, $C{f}_x$ decreases with $\varphi ,{\varphi }_1,{\varphi }_2$ increasing; when $\gamma >-2.50319$, $C{f}_x$ increases with $\varphi ,{\varphi }_1,{\varphi }_2$ increasing. When $\gamma >2.50319$, $C{f}_y$ decreases with $\varphi ,{\varphi }_1,{\varphi }_2$ increasing; when $\gamma <2.50319$, $C{f}_y$ increases with $\varphi ,{\varphi }_1,{\varphi }_2$ increasing. Furthermore, for the axisymmetric case of γ = 0, when the volume fraction $\varphi ={\varphi }_1={\varphi }_2=0$, 5%, 10%, 20%, $C{f}_x=C{f}_y$ = 1.33634, 1.51918, 1.73905, 2.33449. Compared with $\varphi ={\varphi }_1={\varphi }_2=0$, it can be found that the wall shear stress values increase by 13.68%, 30.14%, and 74.69%, respectively. It is shown that the percentage of skin friction coefficient is enhanced by hybrid nanofluid relative to regular fluid. By increasing the fraction of hybrid nanofluids in the fluid, the wall shear stress can be increased. Trends of $C{f}_x$ are opposite for $C{f}_y$. When $\varphi ={\varphi }_1={\varphi }_2$ = 5%, Pr = 6.2, variations of skin friction coefficients $C{f}_x,C{f}_y$ in range of $-6\le \gamma \le 6$ for various σ = 0, 1, 2, 5 are shown in Figure 2b. Moreover, for the axisymmetric case of γ = 0, when volume fraction $\varphi ={\varphi }_1={\varphi }_2$ = 5%, $C{f}_x=C{f}_y$ = 1.51918, 1.87869, 2.18885, 2.93963 for various σ = 0, 1, 2, 5. Compared with $\varphi ={\varphi }_1={\varphi }_2$ = 5%, σ = 0, it can be found that shear stress values increase 23.66%, 40.08%, 93.50%, respectively. It shows that shear stress can be enhanced by increasing the coefficient of permeability in porous media. As σ increases, $C{f}_x$ and $C{f}_y$ increase. When $\varphi$, γ are given certain values, the greater the values of σ, the stronger $C{f}_x$ and $C{f}_y$.

Figure 6a,b illustrates that variations in local Nusselt number $N{u}_x$ against γ in the range of $-6\le \gamma \le 8$ for various $\varphi ={\varphi }_1={\varphi }_2=0$, 5%, 10%, 20%, and σ = 0, 1, 2 when Pr = 6.2. As shown in Figure 3a, when $\varphi ={\varphi }_1={\varphi }_2=0$, 5%, 10%, 20%, σ = 0, and Pr = 6.2, $N{u}_x$ increases as the value of $\varphi ,{\varphi }_1,{\varphi }_2$ increases. When γ < 0, the values of $N{u}_x$ increase steeply with volume fractions $\varphi ,{\varphi }_1,{\varphi }_2$ increasing. When γ > 0, the values of $N{u}_x$ increase and slow down with volume fractions $\varphi ,{\varphi }_1,{\varphi }_2$ increasing. When $\varphi ={\varphi }_1={\varphi }_2$ = 0, 5%, and Pr = 6.2, $N{u}_x$ against γ in the range of $-6\le \gamma \le 8$ for various σ = 0, 1, 2 are shown in Figure 5b. As $\varphi ,{\varphi }_1,{\varphi }_2$ increase, the values of $N{u}_x$ increase. When γ < 0, values of $N{u}_x$ decrease and change clearly with σ increasing. When γ > 0, values of $N{u}_x$ are less affected by the values of σ.

4. Conclusions

Non-axisymmetric Homann stagnant-point flows of the flat plate in porous media containing spherical Cu-Al₂O₃-H₂O nanoparticles are studied using the homotopy analysis method (HAM). The procedure for solving nonlinear coupled ordinary differential equations with variable coefficients involves transforming them into linear ordinary differential equations with constant coefficients. Analytical solutions can be obtained for these equations across the whole domain. A similarity transformation is applied to the governing equations in order to convert them into three coupled nonlinear ordinary differential equations. The analysis involves a discussion of impact of many physical parameters generated in the proposed model. The results have shown that the skin friction coefficients of Cf_x, Cf_y increase with the volume fraction of the hybrid nanofluid and the coefficient of permeability in porous media. For the axisymmetric case of γ = 0, when volume fraction φ = φ₁ = φ₂ = 0, 5%, 10%, 20%, Cf_x = Cf_y = 1.33634, 1.51918, 1.73905, 2.33449. Compared with φ = φ₁ = φ₂ = 0, it can be found that the wall shear stress values increase by 13.68%, 30.14%, 74.69%, respectively. Local Nusselt numbers Nu_x increase as the hybrid nanofluid volume fraction increases. The values of Nu_x decrease, and change clearly with the coefficient of permeability increasing in the range of γ < 0; otherwise, Nu_x are less affected in the range of γ > 0. The shape of nanoparticles has an important effect on the thermal and flow characteristics of nanofluids or hybrid nanofluids. Future works will involve the shape factor of nanoparticles and the mass-based hybrid nanofluid model in addition to the present research.