Electromagnetic Hydrodynamic Convective Flow of Tetra Hybrid Nanofluid in a Porous Medium

1. Introduction

Fluid flow and convective heat transfer are present in a wide range of industrial applications. That is why even today a large amount of research is devoted to these problems, especially to the improvement of heat transfer and management of flow and heat transfer. These requirements are realized in different ways: different external environments such as, for example, electric and magnetic fields, different current spaces such as porous media, different fluids, etc. The idea of improving the thermal conductivity of used fluids is old but its implementation was only started with the publication of the work of Choi and Eastman [1], when fluids in which metal nanoparticles are suspended were practically introduced and called nanofluids (NF).

Among the numerous studies available in the literature that address these problems, only a few of them are mentioned in this introduction. Lima et al. [2] investigated the magnetohydrodynamic (MHD) convective flow of two immiscible fluids in an inclined channel of a porous medium with the effects of buoyancy, Joule and viscous dissipation, wall motion, magnetic field tilt and heat generation/absorption. Raju and Satish [3] investigated the effect of the slip factor on the MHD flow and heat transfer of two incompressible ionized gases inside a horizontal channel with Hall currents. The effects of an inclined magnetic field and thermal radiation on entropy generation during the flow of two immiscible fluids, micropolar and Newtonian, in the porous medium of a horizontal channel were investigated by Yadav et al. [4]. Mixed convective flow of NF in a porous medium in a three-layer vertical channel where the middle layer is a clear fluid and the medium is free was investigated by Umavathi and Sheremet [5], with the influence of an applied magnetic field, heat generation, Brownian motion and heat radiation by Mehta et al. [6] and with a non-Darcy porous medium, an applied magnetic field and four discrete heat sources on the right channel wall Fersadon et al. [7]. Petrović et al. in papers [8] and [9] investigated the EMHD mixed convective flow of two immiscible NF and rigid fluid in porous media in a horizontal and vertical channel, respectively. Subray et al. [10] investigated the EMHD convective flow of NF in a porous medium between two clear fluid layers in an inclined channel with thermal radiation and Hall current. MHD unsteady convective flow of two immiscible NFs in a porous medium in an inclined channel of permeable walls with thermal radiation was investigated by Devi et al. [11].

Kumar and Chamkha [12] investigated Cu/water nanofluid flow in a Darcy-Forchheimer (D-F) porous medium in a convergent/divergent channel with the effects of nanoparticle shape and thermal radiation. Peristaltic flow of NF in a non-Darcy porous medium with thermal radiation, chemical reaction effects and variable viscosity was investigated using the homopathic perturbation method (HPM) by Ouaf and Abouzeid [13]. A review of the literature for various NF models and their solutions-especially accurate, as well as suggestions for further research are given in the work of Sheikh et al. [14]. In a review paper by Nabwey et al. [15] analyzed top papers published between 2018 and 2020 dedicated to the application of NF and heat transfer in porous materials and provided directions for future research. Mixed convective flow of a hybrid nanofluid (HNF) in a vertical channel with slip effects was investigated by Xu and Sun [16]. Zainal et al. [17] investigated the flow and heat transfer characteristics of HNF in the presence of a magnetic field and thermal radiation over a permeable moving surface. Khalil et al. [18] investigated the MHD flow of HNF of variable thermal conductivity and viscosity in a porous medium with heat generation and slip conditions. The effect of the nanoparticle shape factor on the convective flow of immiscible NF and HNF in an inclined channel with a porous medium was investigated by Subray et al. [19]. Comparative analysis of heat and mass transfer of microplate Casson NF (TC4/kerosene oil) and HNF (TC4+NiCr/kerosene oil) with the effects of thermal radiation, magnetic field and porous medium was done by Yahya et al. [20]. Rao and Deka [21] investigated the MHD flow and heat transfer of a water-based HNF induced by a stretchable layer of porous medium with slip boundary conditions and thermal radiation. Mahmood et al. [22] investigated the effect of heat generation/absorption on the MHD flow of three water-based HNFs (THNF) past a stretching/shrinking permeable plate with slip conditions. The effects of radiation, heat generation and nanoparticle size in unsteady MHD mixed convective NF flow in a vertical channel were investigated by Gul et al. [23]. The effects of baffles and chemical reactions on the convective flow of a micropolar fluid in a two-pass vertical channel of a porous medium were investigated by Gitte et al. [24].

MHD flow and heat transfer of a hybrid nanofluid in a boundary layer over an exponentially contracting surface with magnetic field, heat sink/source and thermal radiation effects were investigated numerically by Othman et al. [25]. Kalpana and Saleen [26] investigated the unsteady, laminar flow and heat transfer of a stratified dursty electrically conductive fluid through a porous medium in an infinitely long irregular inclined channel. Agarwal et al. [27] investigated the MHD flow and heat transfer of a micropolar fluid in a permeable channel with thermal radiation using the homotopy perturbation method. The flow and heat transfer of immiscible clear and nanofluids in an inclined channel saturated with a porous medium in the presence of thermal radiation, species diffusion, and viscous and Darcy dissipation effects were investigated by Pavithra et al. [28]. Nikodijević Đorđević et al. [29] investigated the EMHD flow and heat transfer of immiscible nanofluids and hybrid nanofluids in a horizontal channel in which the upper and lower halves of the channel are saturated with different porous media. Parfena et al. [30] considered two-phase Hartmann fluid flow in horizontal and inclined channels where the lower fluid is electrically conductive and the upper complete dielectric in a magnetic field.

Recent developments in mass and heat transport research extensively incorporate TeHNF. Various models are considered including MHD, Soret-Dufour effects, non-Darcy porous media, chemical reactions, heat generation, thermal radiation, Casson fluids, velocity and temperature slip, electroosmosis, blood as a base fluid and others. As can be seen from the works of Amudhini and Poulomi [31], Okasha et al. [32], Mahmood et al. [33], Fatima et al. [34] and Paul and Das [35].

A review of the existing literature indicates an unresolved research gap in this area. To address this gap, the present study focuses on the EMHD mixed convection of TeHNF in a Darcy–Forchheimer medium within a vertical channel, incorporating the effects of thermal radiation. The channel walls are at constant but different temperatures.

2. Mathematical Formulation

This paper analyzes EMHD mixed convection of TeHNF within a vertical channel bounded by two vertical plates and saturated with a porous medium having permeability ${\mathrm{K}}_0$. The porous medium is Darcy-Forchheimer with the Forchheimer inertia factor F and the flow is fully developed. The plates, walls of the channel, are at a distance h from each other and their temperatures are constant. An external homogeneous magnetic field whose induction is B acts perpendicularly to the walls of the channel and a homogeneous electric field whose strength E acts perpendicularly to the plane formed by the directions of the external magnetic field and the primary flow. The pressure gradient in the direction of the primary flow is constant. The coordinate system is chosen so that the X axis is vertical, on the left wall, the Y axis is perpendicular to the channel walls (Figure 1).

When considering this problem, the effects of pressure forces, viscous friction, Darcy’s resistance with Forchheimer’s correction, buoyancy and Lorentz’s, as well as viscous and Darcy-Forchheimer’s dissipation, Joule’s heat and radiant heat are taken into account. For the selected coordinate system, the momentum equation and the energy equation have, respectively, the following forms:

$$-{\displaystyle \frac{\mathrm{d}\mathrm{p}}{\mathrm{d}\mathrm{X}}}+\mathrm{\mu }{\displaystyle \frac{{\mathrm{d}}^2\mathrm{U}}{\mathrm{d}{\mathrm{Y}}^2}}-{\displaystyle \frac{\mathrm{\mu }}{{\mathrm{K}}_0}}\mathrm{U}-\mathrm{\rho }{\displaystyle \frac{\mathrm{F}}{\sqrt{{\mathrm{K}}_0}}}{\mathrm{U}}^2-\mathrm{B}\mathrm{\sigma }\left(\mathrm{E}+\mathrm{B}\mathrm{U}\right)+\mathrm{g}\left(\mathrm{\rho }\mathrm{\beta }\right)\left(\mathrm{T}-{\mathrm{T}}_{{\mathrm{w}}_2}\right)=0,$$

(1)

$$\mathrm{k}{\displaystyle \frac{{\mathrm{d}}^2\mathrm{T}}{\mathrm{d}{\mathrm{Y}}^2}}+\mathrm{\mu }{\left({\displaystyle \frac{\mathrm{d}\mathrm{U}}{\mathrm{d}\mathrm{Y}}}\right)}^2-{\displaystyle \frac{\mathrm{d}{\mathrm{q}}_{\mathrm{r}}}{\mathrm{d}\mathrm{Y}}}+{\displaystyle \frac{\mathrm{\mu }}{{\mathrm{K}}_0}}{\mathrm{U}}^2+\mathrm{\rho }{\displaystyle \frac{\mathrm{F}}{\sqrt{{\mathrm{K}}_0}}}{\mathrm{U}}^3+\mathrm{\sigma }{\left(\mathrm{E}+\mathrm{B}\mathrm{U}\right)}^2=0.$$

(2)

The corresponding boundary conditions are:

$$\left(\mathrm{U},\mathrm{T}\right)\left(0\right)=\left(0,{\mathrm{T}}_{{\mathrm{w}}_2}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(\mathrm{U},\mathrm{T}\right)\left(\mathrm{h}\right)=\left(0,{\mathrm{T}}_{{\mathrm{w}}_1}\right).$$

(3)

The symbols used in equations (1), (2) and the boundary conditions (3) are U, T, g and q_R—fluid velocity, fluid temperature, , gravitational acceleration and radiation heat flux and $\mathrm{\rho },\mathrm{\mu },\mathrm{k},\mathrm{\beta }\hspace{0.17em}\hspace{0.17em}\mathrm{a}\mathrm{n}\mathrm{d}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\sigma }$ are TeHNF properties namely density, dynamic viscosity coefficient, thermal conductivity coefficient, volume thermal expansion coefficient and electrical conductivity coefficient, respectively.

The physical properties of TeHNF are given by the following relations [27]:

$$\mathrm{\mu }={\displaystyle \frac{{\mathrm{\mu }}_{\mathrm{f}}}{\mathrm{m}}},\hspace{0.17em}\hspace{0.17em}\mathrm{k}={\mathrm{\phi }}_1{\mathrm{k}}_{\mathrm{f}},\hspace{0.17em}\hspace{0.17em}\mathrm{\sigma }={\mathrm{\phi }}_2{\mathrm{\sigma }}_{\mathrm{f}},\hspace{0.17em}\hspace{0.17em}\left(\mathrm{\rho }\mathrm{\beta }\right)={\mathrm{\phi }}_3{\left(\mathrm{\rho }\mathrm{\beta }\right)}_{\mathrm{f}},\hspace{0.17em}\hspace{0.17em}\left(\mathrm{\rho }{\mathrm{c}}_{\mathrm{p}}\right)={\mathrm{\phi }}_4{\left(\mathrm{\rho }{\mathrm{c}}_{\mathrm{p}}\right)}_{\mathrm{f}},\hspace{0.17em}\hspace{0.17em}\mathrm{\rho }={\mathrm{\phi }}_5{\mathrm{\rho }}_{\mathrm{f}}$$

(4)

in which, for the sake of brevity, the following labels were used:

$$\begin{array}{l}\mathrm{m}={{\displaystyle \prod _{\mathrm{i}=1}^4\left(1-{\mathrm{\varphi }}_{\mathrm{i}}\right)}}^{2.5},\mathrm{\psi }\left(\mathrm{\varphi },\mathrm{S},\mathrm{R}\right)={\displaystyle \frac{\left(1+2\mathrm{\varphi }\right)\mathrm{S}+2\left(1-\mathrm{\varphi }\right)\mathrm{R}}{\left(1-\mathrm{\varphi }\right)\mathrm{S}+\left(2+\mathrm{\varphi }\right)\mathrm{R}}},\\ {\mathrm{\eta }}_1=\mathrm{\psi }\left({\mathrm{\varphi }}_1,{\mathrm{\sigma }}_1,{\mathrm{\sigma }}_{\mathrm{f}}\right),\hspace{0.17em}\hspace{0.17em}{\mathrm{\sigma }}_{\mathrm{n}\mathrm{f}}={\mathrm{\eta }}_1{\mathrm{\sigma }}_{\mathrm{f}},\hspace{0.17em}\hspace{0.17em}{\mathrm{\eta }}_2=\mathrm{\psi }\left({\mathrm{\varphi }}_2,{\mathrm{\sigma }}_2,{\mathrm{\sigma }}_{\mathrm{n}\mathrm{f}}\right),\\ {\mathrm{\sigma }}_{\mathrm{k}\mathrm{n}\mathrm{f}}={\mathrm{\eta }}_2{\mathrm{\sigma }}_{\mathrm{n}\mathrm{f}},\hspace{0.17em}\hspace{0.17em}{\mathrm{\eta }}_3=\mathrm{\psi }\left({\mathrm{\varphi }}_3,{\mathrm{\sigma }}_3,{\mathrm{\sigma }}_{\mathrm{h}\mathrm{n}\mathrm{f}}\right),\hspace{0.17em}\hspace{0.17em}{\mathrm{\sigma }}_{\mathrm{t}\mathrm{h}\mathrm{n}\mathrm{f}}={\mathrm{\eta }}_3{\mathrm{\sigma }}_{\mathrm{h}\mathrm{n}\mathrm{f}},\\ {\mathrm{\eta }}_4=\mathrm{\psi }\left({\mathrm{\varphi }}_4,{\mathrm{\sigma }}_4,{\mathrm{\sigma }}_{\mathrm{t}\mathrm{h}\mathrm{n}\mathrm{f}}\right),\hspace{0.17em}\hspace{0.17em}{\mathrm{\phi }}_2={\mathrm{\eta }}_1{\mathrm{\eta }}_2{\mathrm{\eta }}_3{\mathrm{\eta }}_4;\hspace{0.17em}\hspace{0.17em}{\mathrm{\phi }}_1={\mathrm{\zeta }}_1{\mathrm{\zeta }}_2{\mathrm{\zeta }}_3{\mathrm{\zeta }}_4;\\ {\mathrm{\zeta }}_1=\mathrm{\psi }\left({\mathrm{\varphi }}_1,{\mathrm{k}}_1,{\mathrm{k}}_{\mathrm{f}}\right),\hspace{0.17em}\hspace{0.17em}{\mathrm{k}}_{\mathrm{n}\mathrm{f}}={\mathrm{\zeta }}_1{\mathrm{k}}_{\mathrm{f}},\hspace{0.17em}\hspace{0.17em}{\mathrm{\zeta }}_2=\mathrm{\psi }\left({\mathrm{\varphi }}_2,{\mathrm{k}}_2,{\mathrm{k}}_{\mathrm{n}\mathrm{f}}\right),\hspace{0.17em}\hspace{0.17em}{\mathrm{k}}_{\mathrm{h}\mathrm{n}\mathrm{f}}={\mathrm{\zeta }}_2{\mathrm{k}}_{\mathrm{n}\mathrm{f}},\\ {\mathrm{\zeta }}_3=\mathrm{\psi }\left({\mathrm{\varphi }}_3,{\mathrm{k}}_3,{\mathrm{k}}_{\mathrm{h}\mathrm{n}\mathrm{f}}\right),\hspace{0.17em}{\mathrm{k}}_{\mathrm{t}\mathrm{h}\mathrm{n}\mathrm{f}}={\mathrm{\zeta }}_3{\mathrm{k}}_{\mathrm{h}\mathrm{n}\mathrm{f}},\hspace{0.17em}\hspace{0.17em}{\mathrm{\zeta }}_4=\mathrm{\psi }\left({\mathrm{\varphi }}_4,{\mathrm{k}}_4,{\mathrm{k}}_{\mathrm{t}\mathrm{h}\mathrm{n}\mathrm{f}}\right);\\ \mathrm{\Omega }\left(\mathrm{H}\right)=\left(1-{\mathrm{\varphi }}_4\right)\left\{\left(1-{\mathrm{\varphi }}_3\right)\left[\left(1-{\mathrm{\varphi }}_2\right)\left(1-{\mathrm{\varphi }}_1+{\mathrm{\varphi }}_1{\displaystyle \frac{{\left(\mathrm{H}\right)}_1}{{\left(\mathrm{H}\right)}_{\mathrm{f}}}}\right)+{\mathrm{\varphi }}_2{\displaystyle \frac{{\left(\mathrm{H}\right)}_2}{{\left(\mathrm{H}\right)}_{\mathrm{f}}}}\right]+{\mathrm{\varphi }}_3{\displaystyle \frac{{\left(\mathrm{H}\right)}_3}{{\left(\mathrm{H}\right)}_{\mathrm{f}}}}\right\}+{\mathrm{\varphi }}_4{\displaystyle \frac{{\left(\mathrm{H}\right)}_4}{{\left(\mathrm{H}\right)}_{\mathrm{f}}}},\\ {\mathrm{\phi }}_3=\mathrm{\Omega }\left(\left(\mathrm{\rho }{\mathrm{\beta }}_{\mathrm{t}}\right)\right),\hspace{0.17em}\hspace{0.17em}{\mathrm{\phi }}_3=\mathrm{\Omega }\left(\left(\mathrm{\rho }{\mathrm{c}}_{\mathrm{p}}\right)\right),\hspace{0.17em}\hspace{0.17em}{\mathrm{\phi }}_5=\mathrm{\Omega }\left(\mathrm{\rho }\right)\end{array}$$

(5)

where the subscripts 1, 2, 3 and 4 in the physical properties refer to the materials of nanoparticles 1, 2, 3 and 4 and their volume fractions and the subscript f to the properties of the base fluid.

Rasseland transformation [4] is used for radiative heat flux:

$${\mathrm{q}}_{\mathrm{r}}=-{\displaystyle \frac{4{\mathrm{\sigma }}^{*}}{3{\mathrm{k}}^{*}}}{\displaystyle \frac{\partial {\mathrm{T}}^4}{\partial \mathrm{Y}}}$$

(6)

in which σ* is the Stefan-Boltzmann constant and k* is the mean value of the absorption coefficient of TeHNF. Expanding T₄ in order and keeping to two terms, is obtained:

$${\displaystyle \frac{\mathrm{d}{\mathrm{q}}_{\mathrm{r}}}{\mathrm{d}\mathrm{Y}}}=-{\displaystyle \frac{16{\mathrm{\sigma }}^{*}{\mathrm{T}}_{{\mathrm{w}}_2}^3}{3{\mathrm{k}}^{*}}}{\displaystyle \frac{{\mathrm{d}}^2\mathrm{T}}{\mathrm{d}{\mathrm{Y}}^2}}.$$

(7)

For further study of the considered problem, it is more convenient to transform equations (1) and (2), as well as boundary conditions (3), into dimensionless quantities. For this purpose, the following dimmensionless quantities have been introduced:

$$\mathrm{y}={\displaystyle \frac{\mathrm{Y}}{\mathrm{h}}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{u}={\displaystyle \frac{\mathrm{U}}{{\mathrm{U}}_0}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\theta }={\displaystyle \frac{\mathrm{T}-{\mathrm{T}}_{{\mathrm{w}}_2}}{{\mathrm{T}}_{{\mathrm{w}}_1}-{\mathrm{T}}_{{\mathrm{w}}_2}}},$$

(8)

where Uo is the characteristic velocity, which will be chosen in the following text.

By replacing physical properties (5), relations (6) and quantities (8) in equations (1) and (2) as well as boundary conditions (3), the equations are obtained, respectively:

$${\displaystyle \frac{{\mathrm{d}}^2\mathrm{u}}{\mathrm{d}{\mathrm{y}}^2}}-{\mathrm{\omega }}^2\mathrm{u}-{\mathrm{R}}_1+\mathrm{a}\mathrm{M}\mathrm{\theta }-\mathrm{f}{\mathrm{u}}^2=0,$$

(9)

$${\displaystyle \frac{{\mathrm{d}}^2\mathrm{\theta }}{\mathrm{d}{\mathrm{y}}^2}}+\mathrm{n}\left[{\left({\displaystyle \frac{\mathrm{d}\mathrm{u}}{\mathrm{d}\mathrm{y}}}\right)}^2+\mathrm{f}{\mathrm{u}}^3+{\mathrm{\omega }}^2{\mathrm{u}}^2+2\mathrm{R}\mathrm{u}+{\mathrm{R}}_2\right]=0$$

(10)

and boundary conditions:

$$(\mathrm{u},\mathrm{\theta })(0)=(0,0),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(\mathrm{u},\mathrm{\theta })(1)=(0,1).$$

(11)

The following notations were used in the last equations and boundary conditions:

$$\begin{array}{l}\mathrm{\Lambda }={\displaystyle \frac{{\mathrm{h}}^2}{{\mathrm{K}}_0}},\hspace{0.17em}\hspace{0.17em}\mathrm{Re}={\displaystyle \frac{\mathrm{h}{\mathrm{U}}_0}{{\mathrm{\mu }}_{\mathrm{f}}}}{\mathrm{\rho }}_{\mathrm{f}},\hspace{0.17em}\hspace{0.17em}{\mathrm{\Lambda }}_{\mathrm{f}}={\displaystyle \frac{\mathrm{h}\mathrm{F}}{\sqrt{{\mathrm{K}}_{\mathrm{o}}}}},\hspace{0.17em}\hspace{0.17em}\mathrm{H}\mathrm{a}=\mathrm{h}\mathrm{B}\sqrt{{\displaystyle \frac{{\mathrm{\sigma }}_{\mathrm{f}}}{{\mathrm{\mu }}_{\mathrm{f}}}},}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{K}={\displaystyle \frac{\mathrm{E}}{\mathrm{B}{\mathrm{U}}_0}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{n}={\displaystyle \frac{\mathrm{B}\mathrm{r}}{\mathrm{m}\mathrm{b}}},\\ \mathrm{M}=\mathrm{g}{\displaystyle \frac{{\left(\mathrm{\rho }\mathrm{\beta }\right)}_{\mathrm{f}}{\mathrm{h}}^2\left({\mathrm{T}}_{{\mathrm{w}}_1}-{\mathrm{T}}_{{\mathrm{w}}_2}\right)}{{\mathrm{\mu }}_{\mathrm{f}}{\mathrm{U}}_0}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{P}=-{\displaystyle \frac{\partial \mathrm{p}}{\partial \mathrm{X}}}{\displaystyle \frac{{\mathrm{h}}^2}{{\mathrm{\mu }}_{\mathrm{f}}{\mathrm{U}}_0}},\hspace{0.17em}\hspace{0.17em}{\mathrm{\omega }}^2=\mathrm{\Lambda }+\mathrm{m}{\mathrm{\phi }}_2\mathrm{H}{\mathrm{a}}^2,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{a}=\mathrm{m}{\mathrm{\phi }}_3,\\ \mathrm{f}=\mathrm{m}{\mathrm{\phi }}_5\mathrm{Re}{\mathrm{\Lambda }}_{\mathrm{f}},\hspace{0.17em}\hspace{0.17em}\mathrm{R}=\mathrm{m}{\mathrm{\phi }}_2\mathrm{K}\mathrm{H}{\mathrm{a}}^2,\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}_1=\mathrm{R}-\mathrm{m}\mathrm{P},\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}^{*}={\displaystyle \frac{16{\mathrm{\sigma }}^{*}{\mathrm{T}}_{{\mathrm{w}}_2}^3}{3{\mathrm{k}}^{*}{\mathrm{k}}_{\mathrm{f}}}},\hspace{0.17em}\hspace{0.17em}\mathrm{b}={\mathrm{\phi }}_1+{\mathrm{R}}^{*},\\ \mathrm{Pr}={\displaystyle \frac{{\mathrm{\mu }}_{\mathrm{f}}{\mathrm{c}}_{\mathrm{p}\mathrm{f}}}{{\mathrm{k}}_{\mathrm{f}}}},\hspace{0.17em}\hspace{0.17em}\mathrm{E}\mathrm{c}={\displaystyle \frac{{\mathrm{U}}_0^2}{{\mathrm{c}}_{\mathrm{p}\mathrm{f}}\left({\mathrm{T}}_{{\mathrm{w}}_1}-{\mathrm{T}}_{{\mathrm{w}}_2}\right)}},\hspace{0.17em}\hspace{0.17em}\mathrm{B}\mathrm{r}=\mathrm{Pr}\mathrm{E}\mathrm{c},\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}_2=\mathrm{K}\mathrm{R}\end{array}$$

(12)

where: Λ-porosity factor, K-factor of external electric load factor and Reynolds (Re), Hartmann (Ha), Prandtl (Pr), Eckert (Ec) and Brinkman (Br) numbers.

Equations (9) and (10), as well as boundary conditions (11) are dimensionless and represent a mathematical model of the mixed convective heat transfer problem considered here.

3. Solutions

Here, to determine the solution to the described problem, i.e., to solve equations (9) and (10) with boundary conditions (11) he uses the Homotopy perturbation method (HPM) authored by J. H. He [36]. Following this method, homotopy equations corresponding to equations (9) and (10) are formed, respectively.

$${\displaystyle \frac{{\mathrm{d}}^2\mathrm{u}}{\mathrm{d}{\mathrm{y}}^2}}-{\mathrm{\omega }}^2\mathrm{u}-{\mathrm{R}}_1+\mathrm{\lambda }\left(\mathrm{a}\mathrm{M}\mathrm{\theta }-\mathrm{f}{\mathrm{u}}^2\right)=0,$$

(13)

$${\displaystyle \frac{{\mathrm{d}}^2\mathrm{\theta }}{\mathrm{d}{\mathrm{y}}^2}}+\mathrm{S}+\mathrm{\lambda }\mathrm{n}\left[{\left({\displaystyle \frac{\mathrm{d}\mathrm{u}}{\mathrm{d}\mathrm{y}}}\right)}^2+\mathrm{f}{\mathrm{u}}^3+{\mathrm{\omega }}^2{\mathrm{u}}^2+2\mathrm{R}\mathrm{u}\right]=0$$

(14)

where S=nR₂ and $\mathrm{\lambda }\in \left[0,1\right]$ the homotopy parameter. The solution of the last equations is assumed in the form:

$$(\mathrm{u},\mathrm{\theta })(\mathrm{y},\mathrm{\lambda })=({\mathrm{u}}_0,{\mathrm{\theta }}_0)(\mathrm{y})+\mathrm{\lambda }({\mathrm{u}}_1,{\mathrm{\theta }}_1)(\mathrm{y})+\mathrm{...}$$

(15)

Substituting relations (15) in equations (13) and (14) yields zero-order equations:

$${\displaystyle \frac{{\mathrm{d}}^2{\mathrm{u}}_0}{\mathrm{d}{\mathrm{y}}^2}}-{\mathrm{\omega }}^2{\mathrm{u}}_0={\mathrm{R}}_1,$$

(16)

$${\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\theta }}_0}{\mathrm{d}{\mathrm{y}}^2}}=-\mathrm{S},$$

(17)

and first-order equations:

$${\displaystyle \frac{{\mathrm{d}}^2{\mathrm{u}}_1}{\mathrm{d}{\mathrm{y}}^2}}-{\mathrm{\omega }}^2{\mathrm{u}}_1+\mathrm{a}\mathrm{M}{\mathrm{\theta }}_0-\mathrm{f}{\mathrm{u}}_0^2=0,$$

(18)

$${\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\theta }}_1}{\mathrm{d}{\mathrm{y}}^2}}+\mathrm{n}\left[{\left({\displaystyle \frac{\mathrm{d}{\mathrm{u}}_0}{\mathrm{d}\mathrm{y}}}\right)}^2+\mathrm{f}{\mathrm{u}}_0^3+{\mathrm{\omega }}^2{\mathrm{u}}_0^2+2\mathrm{R}{\mathrm{u}}_0\right]=0.$$

(19)

The corresponding zero and first order boundary conditions are:

$$({\mathrm{u}}_0,{\mathrm{\theta }}_0)(0)=(0,0),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}({\mathrm{u}}_0,{\mathrm{\theta }}_0)(1)=(0,1),$$

(20)

$$({\mathrm{u}}_1,{\mathrm{\theta }}_1)(0)=(0,0),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}({\mathrm{u}}_1,{\mathrm{\theta }}_1)(1)=(0,0).$$

(21)

respectively.

The solution of the zero-order equations that satisfies the appropriate boundary conditions (20) is given by the following relations:

$${\mathrm{u}}_0\left(\mathrm{y}\right)={\mathrm{C}}_1\exp \left(\mathrm{\omega }\mathrm{y}\right)+{\mathrm{C}}_2\exp \left(-\mathrm{\omega }\mathrm{y}\right)-\mathrm{D},$$

(22)

$${\mathrm{\theta }}_0\left(\mathrm{y}\right)=-{\displaystyle \frac{\mathrm{S}}{2}}{\mathrm{y}}^2+{\mathrm{D}}_1\mathrm{y},$$

(23)

in which the following tags were used:

$$\mathrm{D}={\displaystyle \frac{{\mathrm{R}}_1}{{\mathrm{\omega }}^2}},\hspace{0.17em}\hspace{0.17em}{\mathrm{C}}_1=\mathrm{D}{\displaystyle \frac{1-\exp \left(-\mathrm{\omega }\right)}{\exp \left(\mathrm{\omega }\right)-\exp \left(-\mathrm{\omega }\right)}},\hspace{0.17em}\hspace{0.17em}{\mathrm{C}}_2=\mathrm{D}-{\mathrm{C}}_1,\hspace{0.17em}\hspace{0.17em}{\mathrm{D}}_1=1+{\displaystyle \frac{\mathrm{S}}{2}}.$$

(24)

Substituting relations (22) and (23) in the first-order equations and solving them shows that their solution, which satisfies the boundary conditions (21), is given by the relations:

$$\begin{array}{l}{\mathrm{u}}_1\left(\mathrm{y}\right)={\mathrm{C}}_3\exp \left(\mathrm{\omega }\mathrm{y}\right)+{\mathrm{C}}_4\exp \left(-\mathrm{\omega }\mathrm{y}\right)+{\mathrm{R}}_3\exp \left(2\mathrm{\omega }\mathrm{y}\right)+{\mathrm{R}}_4\exp \left(-2\mathrm{\omega }\mathrm{y}\right)-\\ {\mathrm{R}}_5\mathrm{y}\exp \left(\mathrm{\omega }\mathrm{y}\right)+{\mathrm{R}}_6\mathrm{y}\exp \left(-\mathrm{\omega }\mathrm{y}\right)+{\mathrm{a}}_1{\mathrm{y}}^2+{\mathrm{a}}_2\mathrm{y}+{\mathrm{a}}_3,\end{array}$$

(25)

$$\begin{array}{l}{\mathrm{\theta }}_1\left(\mathrm{y}\right)=-\mathrm{n}\left[{\mathrm{R}}_9\exp \left(3\mathrm{\omega }\mathrm{y}\right)+{\mathrm{R}}_{10}\exp \left(-3\mathrm{\omega }\mathrm{y}\right)+{\mathrm{R}}_{11}\exp \left(2\mathrm{\omega }\mathrm{y}\right)+\right.\\ {\mathrm{R}}_{12}\exp \left(-2\mathrm{\omega }\mathrm{y}\right)+\left.{\mathrm{R}}_{13}\exp \left(\mathrm{\omega }\mathrm{y}\right)+{\mathrm{R}}_{14}\exp \left(-\mathrm{\omega }\mathrm{y}\right)+{\mathrm{R}}_{15}{\mathrm{y}}^2+{\mathrm{D}}_3\mathrm{y}+{\mathrm{D}}_4\right]\end{array}$$

(26)

in which the following tags were used:

$$\begin{array}{l}{\mathrm{a}}_1=-{\displaystyle \frac{\mathrm{a}\mathrm{M}\mathrm{S}}{2{\mathrm{\omega }}^2}},\hspace{0.17em}\hspace{0.17em}{\mathrm{a}}_2=-{\displaystyle \frac{\mathrm{a}\mathrm{M}{\mathrm{D}}_1}{{\mathrm{\omega }}^2}},\hspace{0.17em}\hspace{0.17em}{\mathrm{a}}_3={\displaystyle \frac{1}{{\mathrm{\omega }}^2}}\left[2{\mathrm{a}}_1-\mathrm{f}\left({\mathrm{D}}^2+2{\mathrm{C}}_1{\mathrm{C}}_2\right)\right],\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}_9={\displaystyle \frac{\mathrm{f}{\mathrm{C}}_1^3}{9{\mathrm{\omega }}^2}},\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}_{10}={\displaystyle \frac{\mathrm{f}{\mathrm{C}}_2^3}{9{\mathrm{\omega }}^2}},\\ {\mathrm{R}}_3=\mathrm{f}{\displaystyle \frac{{\mathrm{C}}_1^2}{3{\mathrm{\omega }}^2}},\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}_4=\mathrm{f}{\displaystyle \frac{{\mathrm{C}}_2^2}{3{\mathrm{\omega }}^2}},\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}_5=\mathrm{f}{\displaystyle \frac{\mathrm{D}{\mathrm{C}}_1}{\mathrm{\omega }}},\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}_6=\mathrm{f}{\displaystyle \frac{\mathrm{D}{\mathrm{C}}_2}{\mathrm{\omega }}},\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}_7={\mathrm{R}}_3+{\mathrm{R}}_4+{\mathrm{a}}_3,\\ {\mathrm{R}}_8={\mathrm{R}}_3\exp \left(2\mathrm{\omega }\right)+{\mathrm{R}}_4\exp \left(-2\mathrm{\omega }\right)-{\mathrm{R}}_5\exp \left(\mathrm{\omega }\right)+{\mathrm{R}}_6\exp \left(-\mathrm{\omega }\right)+{\mathrm{a}}_1+{\mathrm{a}}_2+{\mathrm{a}}_3,\\ {\mathrm{C}}_3={\displaystyle \frac{{\mathrm{R}}_7\exp \left(-\mathrm{\omega }\right)-{\mathrm{R}}_8}{\exp \left(\mathrm{\omega }\right)-\exp \left(-\mathrm{\omega }\right)}},\hspace{0.17em}\hspace{0.17em}{\mathrm{C}}_4=-{\mathrm{C}}_3-{\mathrm{R}}_7,\hspace{0.17em}\hspace{0.17em}\mathrm{d}={\displaystyle \frac{1}{4{\mathrm{\omega }}^2}}\left(2{\mathrm{\omega }}^2-3\mathrm{f}\mathrm{D}\right),\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}_{11}=\mathrm{d}{\mathrm{C}}_1^2,\hspace{0.17em}\\ {\mathrm{R}}_{12}=\mathrm{d}{\mathrm{C}}_2^2,\hspace{0.17em}\hspace{0.17em}{\mathrm{d}}_1={\displaystyle \frac{1}{{\mathrm{\omega }}^2}}\left[3\mathrm{f}\left({\mathrm{C}}_1{\mathrm{C}}_2+{\mathrm{D}}^2\right)+2\left(\mathrm{R}-{\mathrm{\omega }}^2\mathrm{D}\right)\right],\\ {\mathrm{R}}_{13}={\mathrm{d}}_1{\mathrm{C}}_1,{\mathrm{R}}_{14}={\mathrm{d}}_1{\mathrm{C}}_2,\\ {\mathrm{R}}_{15}={\displaystyle \frac{1}{2}}\left[\left({\mathrm{\omega }}^2\mathrm{D}-2\mathrm{R}-\mathrm{f}{\mathrm{D}}^2\right)\mathrm{D}-2{\mathrm{C}}_1{\mathrm{C}}_2\left({\mathrm{\omega }}^2+3\mathrm{f}{\mathrm{D}}_1\right)\right],\\ {\mathrm{R}}_{16}={\mathrm{R}}_9+{\mathrm{R}}_{10}+{\mathrm{R}}_{11}+{\mathrm{R}}_{12}+{\mathrm{R}}_{13}+{\mathrm{R}}_{14},\\ {\mathrm{D}}_3={\mathrm{R}}_{16}-{\mathrm{R}}_{15}-{\mathrm{R}}_9\exp \left(3\mathrm{\omega }\right)-{\mathrm{R}}_{10}\exp \left(-3\mathrm{\omega }\right)-{\mathrm{R}}_{11}\exp \left(2\mathrm{\omega }\right)-{\mathrm{R}}_{12}\exp \left(-2\mathrm{\omega }\right)\\ -{\mathrm{R}}_{13}\exp \left(\mathrm{\omega }\right)-{\mathrm{R}}_{14}\exp \left(-\mathrm{\omega }\right),\hspace{0.17em}\hspace{0.17em}{\mathrm{D}}_4=-{\mathrm{R}}_9-{\mathrm{R}}_{10}-{\mathrm{R}}_{11}-{\mathrm{R}}_{12}-{\mathrm{R}}_{13}-{\mathrm{R}}_{14}.\end{array}$$

(27)

The solution of equations (9) and (10) is obtained from relations (15) when $\mathrm{\lambda }\to 1$ it is also given by the relations:

$$(\mathrm{u},\mathrm{\theta })(\mathrm{y})=({\mathrm{u}}_0,{\mathrm{\theta }}_0)(\mathrm{y})+({\mathrm{u}}_1,{\mathrm{\theta }}_1)(\mathrm{y})+\mathrm{...}$$

(28)

where u₀, θ₀ and u₁,θ₁ are already determined quantities.

Relations (28) represent the velocity and temperature distribution of the TeHNF problem described here.

For the analysis of this problem, the shear stresses on the channel walls, as well as the local Nusselt numbers on these walls, are of interest. The shear stresses on the left and right channel walls are given by the relations:

$${\mathrm{\tau }}_{1,2}={\displaystyle \frac{1}{\mathrm{m}}}{\left.{\displaystyle \frac{\mathrm{d}\mathrm{u}}{\mathrm{d}\mathrm{y}}}\right|}_{\mathrm{y}=0,1},$$

(29)

respectively. The Nusselt numbers on these walls are given by the relations:

$$\mathrm{N}{\mathrm{u}}_{1,2}={\mathrm{\phi }}_1{\left.{\displaystyle \frac{\mathrm{d}\mathrm{\theta }}{\mathrm{d}\mathrm{y}}}\right|}_{\mathrm{y}=0,1},$$

(30)

respectively.

The characteristic velocity U₀ has not yet been selected, and this can be done in different ways. Here it will be chosen so that the quantity M becomes Grashof’s number, i.e., from the equation:

$$\mathrm{g}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{f}}{\mathrm{\beta }}_{\mathrm{f}}{\mathrm{h}}^2\left({\mathrm{T}}_{{\mathrm{w}}_1}-{\mathrm{T}}_{{\mathrm{w}}_2}\right)}{{\mathrm{U}}_0{\mathrm{\mu }}_{\mathrm{f}}}}=\mathrm{g}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{f}}^2{\mathrm{\beta }}_{\mathrm{f}}{\mathrm{h}}^3\left({\mathrm{T}}_{{\mathrm{w}}_1}-{\mathrm{T}}_{{\mathrm{w}}_2}\right)}{{\mathrm{\mu }}_{\mathrm{f}}^2}},$$

(31)

from which it is obtained that:

$${\mathrm{U}}_0={\displaystyle \frac{{\mathrm{\mu }}_{\mathrm{f}}}{\mathrm{h}{\mathrm{\rho }}_{\mathrm{f}}}}.$$

(32)

4. Model Verification

In order to verify the flow and heat transfer model considered here, the results obtained in this work are compared with the results obtained by Umavathi and Malashetty [37]. For the case where ${\mathrm{\varphi }}_1={\mathrm{\varphi }}_2={\mathrm{\varphi }}_3={\mathrm{\varphi }}_4=0,\hspace{0.17em}\mathrm{H}\mathrm{a}=2,\hspace{0.17em}\mathrm{\Lambda }=0,\hspace{0.17em}\mathrm{F}=0,\hspace{0.17em}\mathrm{K}=0,\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}^{*}=0\hspace{0.17em}$the problem considered here is reduced to the problem considered by Umavathi and Malashetty [37]. The results obtained here for the fluid velocity distribution and the results obtained in [37] are presented in Figure 2. This figure shows that results match quite well, i.e., that the differences between them are minimal, do not exceed 2%, which is quite satisfactory.

5. Results and Analysis

In the previous chapter, the distributions of velocity and temperature, as well as the shear stresses and the local Nusselt numbers on the channel walls, were determined as functions of the introduced physical parameters. Part of the obtained results for the physical properties of the nanoparticles and the base fluid is presented in Table 1, while the results corresponding to different values of the introduced parameters are given in this chapter in the form of graphs and Table 2. Specifically, the results are presented for the case in which the base fluid is water and the nanoparticles are aluminum oxide (Al₂O₃), titanium dioxide (TiO₂), magnesium oxide (MgO) and magnetite (Fe₃O₄). The values of the physical parameters used in the present analysis are: $\mathrm{P}=25,\hspace{0.17em}\hspace{0.17em}\mathrm{F}=0.01,\hspace{0.17em}\hspace{0.17em}\mathrm{H}\mathrm{a}=2,\hspace{0.17em}\hspace{0.17em}$$\mathrm{K}=-1,\hspace{0.17em}\hspace{0.17em}\mathrm{M}=0.5,\hspace{0.17em}\hspace{0.17em}{\mathrm{R}}^{*}=0.6,\hspace{0.17em}\hspace{0.17em}\mathrm{B}\mathrm{r}=0.5,\hspace{0.17em}\hspace{0.17em}\mathrm{\Lambda }=5\hspace{0.17em}\hspace{0.17em}\mathrm{a}\mathrm{n}\mathrm{d}\hspace{0.17em}\hspace{0.17em}{\mathrm{\varphi }}_1={\mathrm{\varphi }}_2={\mathrm{\varphi }}_3={\mathrm{\varphi }}_4=0.01$, except when these parameters are varied to investigate their individual effects.

Figure 3 and Figure 4 show the distributions of fluid velocity and temperature in the channel for different values of the Hartmann number (Ha), respectively. These figures indicate that higher values of Ha correspond to lower fluid velocities and higher temperatures. This behavior can be explained by the fact that increasing the Hartmann number strengthens the Lorentz force, which in this case acts as a resistive force, thereby reducing the fluid velocity and its gradient in the vicinity of the channel walls. The enhanced Lorentz force also increases Joule dissipation, which leads to a rise in the fluid temperature within the channel.

The influence of the porosity parameter (Λ) on the velocity and temperature distributions of the fluid in the channel is illustrated in Figure 5 and Figure 6, respectively. These figures show that higher values of the porosity parameter correspond to lower fluid velocities and higher temperatures within the channel. Larger values of Λ represent porous media with lower permeability, i.e., media with a denser porous matrix. As the fluid flows through a denser matrix, the Darcy and Forchheimer resistances to the fluid motion increase, resulting in a reduction in fluid velocity. During the passage of the fluid through the porous medium, the energy expended in overcoming this resistance is converted into internal thermal energy, which leads to an increase in the fluid temperature.

Figure 7 and Figure 8 show the distributions of fluid velocity and temperature in the channel for different values of the Forchheimer parameter (F), respectively. These figures indicate that increasing the value of F reduces the fluid velocity and causes a slight change in the fluid temperature within the channel. This behavior can be explained by the fact that an increase in F enhances the Forchheimer (inertial) drag force of the porous medium, which leads to a reduction in the fluid flow velocity. This reduction in velocity has a dominant effect on the energy balance: it decreases heat generation due to Forchheimer, Darcy, viscous and Joule dissipation, thereby leading to a decrease in the fluid temperature in the channel.

The velocity distributions of the fluid for different values of the Grashof number (M) are shown in Figure 9. It can be observed that higher values of M correspond to higher fluid velocities. This result can be explained by the fact that larger values of M represent stronger buoyancy forces, which enhance the fluid motion. In the present case, the parameter M has a primarily dynamical rather than thermal influence; therefore, the temperature field remains essentially unchanged and is not presented here.

The effects of the parameter P on the velocity and temperature distributions of the fluid in the channel are illustrated in Figure 10 and Figure 11, respectively. These figures show that increasing the value of P enhances the fluid velocity and leads to a rise in the fluid temperature within the channel. An increase in P strengthens the driving force due to the larger pressure gradient, which in turn increases the flow velocity. The higher velocity intensifies the dissipation of flow energy in the channel, resulting in an increase in the fluid temperature.

Figure 12 presents the temperature distributions of the fluid in the channel for different values of the radiation parameter (R*). It can be observed that increasing R* reduces the fluid temperature due to the enhanced heat transfer by thermal radiation. In this case, the dominant effects on the fluid flow are the pressure gradient and the magnetic field rather than the buoyancy effect, therefore, the velocity field remains highly stable and is not presented here.

The effects of the external electric loading parameter (K) on the velocity and temperature distributions are shown in Figure 13 and Figure 14, respectively. The sign of this parameter depends on the direction of the applied external electric field. For negative values of K, the electric-field-induced force acting on the fluid is a driving force, whereas for positive values it acts as a resistive force opposing the fluid motion. As shown in Figure 13, when K = − 0.1 the fluid velocity increases, while for K = 0.1 the velocity decreases compared to the case K = 0, which is also reflected in the corresponding temperature distributions in Figure 14.

Figure 15 and Figure 16 show the velocity and temperature distributions of the fluid in the channel for different values of the Brinkman number (Br), respectively. Higher values of Br correspond to larger amounts of heat generated by viscous dissipation in the fluid and consequently to higher temperatures, as observed in Figure 16. The increased temperature enhances the buoyancy force, which, despite the presence of resistive forces, can increase the fluid velocity in the channel, as illustrated in Figure 15.

Figure 17 and Figure 18 present the velocity and temperature distributions for a clear fluid, a nanofluid, a hybrid nanofluid, a trihybrid nanofluid and a tetra-hybrid nanofluid, respectively. The properties of these fluids—density, dynamic viscosity, thermal conductivity and electrical conductivity become progressively more pronounced when moving from the clear fluid to the tetra-hybrid nanofluid. From these figures, it is observed that both the velocity and the temperature of the fluid in the channel decrease as one moves from the clear fluid to the tetra-hybrid nanofluid.

These trends can be physically explained as follows. An increase in fluid viscosity reduces its kinematic mobility. An increase in density raises the inertia, which, in combination with the Forchheimer force, further slows down the fluid motion. An increase in electrical conductivity strengthens the Lorentz force, thereby further slowing down the flow. In addition, Darcy resistance also contributes to flow deceleration. All these effects lead to a reduction in the fluid velocity when moving from the clear fluid to the tetra-hybrid nanofluid.

The decrease in velocity reduces viscous, Darcy and Forchheimer dissipation. Moreover, the increase in thermal conductivity enhances heat transfer toward the channel walls, which are maintained at constant temperatures. Thermal radiation further promotes heat removal. Although higher electrical conductivity increases Joule dissipation, in the present case its effect is outweighed by the aforementioned mechanisms. As a result, less heat remains in the fluid when moving from the clear fluid to the tetra-hybrid nanofluid, leading to a decrease in temperature, as observed in Figure 18.

The values of the wall shear stresses and the local Nusselt numbers at the channel walls are listed in Table 2.

Table 2 gives numerical values of shear stress and Nusselt numbers on channel walls. From this table, it can be seen that the increase in the values of the parameters Ha, Λ, F and R* leads to a decrease in the value of the shear stresses on the channel walls. An increase in the values of the parameters M, P and Br leads to an increase in the values of the shear stresses on the channel walls. An increase in the value of the parameters F and R* leads to a decrease in the value of the Nusselt numbers in the channels, and the change in the value of M in the walls does not affect their values. An increase in the parameters Ha, Λ, P, R* and Br leads to an increase in the value of Nusselt numbers.

6. Conclusions

In this paper, electromagnetic hydrodynamic flow and heat transfer of tetra hybrid nanofluid in Darcy-Forchheimer porous medium in a vertical channel were investigated. Darcy, Forchheimer, Joule and radiant heat dissipations are included. The resulting dimensionless equations were solved using the homotopy perturbation method. Analytical distributions of fluid velocity and temperature, as well as shear stresses and Nusselt numbers on channel walls, were determined as a function of the introduced physical parameters significant for this problem. Part of the specific numerical results for velocity and temperature distributions is given in the form of graphs, and for shear stresses and Nusselt numbers on channel walls in a table.

Based on the analysis of the obtained results, certain conclusions were drawn, some of which are given in the following text.

An increase in the value of the Ha number, as well as an increase in the value of the porosity factor, cause a decrease in the fluid velocity and shear stress on the channel walls, and an increase in the fluid temperature and an increase in the value of Nusselt numbers on the channel walls.

Higher values of the Forchheimer factor, as well as higher values of heat radiation, correspond to lower velocities, lower temperatures, lower values of shear stresses and Nusselt numbers on the channel walls.

An increase in the value of the pressure gradient as well as an increase in the value of the Brinkman number cause the growth of all the previously mentioned analyzed quantities $\left(\mathrm{u},\mathrm{T},{\mathrm{\tau }}_1,{\mathrm{\tau }}_2,{\mathrm{N}}_{{\mathrm{u}}_1},{\mathrm{N}}_{{\mathrm{u}}_2}\right).$An increase in the Grashof number causes an increase in the velocity and shear stress on the channel walls, and has no effect on the fluid temperature and Nusselt numbers on the channel walls.

The direction and strength of the external electric field, i.e., the sign and value of the external electric load factor influence the increase/decrease of the analyzed quantities in relation to their values when this field is absent.

Tetra hybrid nanofluid shows advantages over simpler fluids with fewer or no nanoparticle species.