Computation of non-Newtonian quadratic convection in electro-magneto-hydrodynamic (EMHD) duct flow with temperature-dependent viscosity

Abstract

Motivated by novel developments in smart non-Newtonian thermal duct systems, a theoretical study has been presented in this article for electro-magneto-hydrodynamic (EMHD) buoyancy-driven flow of a fourth-grade viscoelastic fluid in a vertical duct with quadratic convection. The viscosity of the fourth-grade fluid model is assumed to be temperature-dependent, and the Reynolds exponential model is deployed. Viscous heating and Joule dissipation effects are included. The duct comprises a pair of parallel electrically insulated vertical flat plates located a finite distance apart. Via suitable scaling transformations, a nonlinear boundary value problem is derived for the momentum and heat transport. A homotopy perturbation method (HPM) solution is obtained coded in Mathematica symbolic software. There is a considerable enhancement in wall skin friction with an increment in fourth-grade fluid parameter, Brinkman number, electrical field parameter, thermal buoyancy parameter, and quadratic thermal convection parameter. However, skin friction is strongly reduced with a rise in variable viscosity parameter, Hartmann (magnetic) number, and electromagnetic heat generation to conduction ratio. Nusselt number magnitudes are elevated with increase in variable viscosity parameter, thermal buoyancy parameter, and quadratic thermal convection parameter, whereas they are significantly suppressed with increment in fourth-grade fluid parameter, Brinkman number, and Hartmann magnetic number.

Introduction

Magnetohydrodynamics (MHD) has emerged as a major discipline in 21st century engineering sciences, tracing its roots to astrophysics. ¹ It involves the interaction between magnetic fields and either viscous or inviscid electrically conducting fluids. In modern energy systems a major application is the heat transfer and flow control in nuclear ducts and other energy systems (solar) where excessive temperatures can arise. MHD ducts flows feature in lithium blankets, ² turbulence suppression, ³ magnetic sensors, ⁴ ferromagnetic solar absorber collectors, ⁵ liquid metal flow couplers, ⁶ coolant deployment in sodium fast reactors (SFRs), ⁷ and pressure drop regulation. ⁸ In mathematical models of MHD duct flows frequently dimensionless parameters arise which include the Hartmann number (ratio of magnetic body force to viscous force), Stuart magnetic interaction number (ratio of magnetic force to inertial force), and magnetic Reynolds number (when induction effects are significant). Electrohydrodynamics (EHD) ⁹ is a sub-discipline of the more general continuum electromechanics and also known as electro-kinetics. It focuses on the interaction between electrical fields and electrically charged liquids. While in MHD usually transverse magnetic fields are utilized to generate a perpendicular Lorentz electromagnetic body force, in EHD the electrical field is often implemented along the axis of the conduit or duct. Many exciting applications for EHD duct flows (e.g. parallel plate geometries) have emerged in recent decades including atomizers, ¹⁰ corona-based duct flows, ¹¹ electrode gas pump channels, ¹² electrostatic precipitators, ¹³ dielectric barriers for electro-fluid injection, ¹⁴ EHD heat exchangers, ¹⁵ coolants for circuits, ¹⁶ thermal enhancement, ¹⁷ and bio-inspired microfluidic pumping systems. ¹⁸ As with MHD, some popular electrohydrodynamic dimensionless numbers arise in numerical and theoretical studies of EHD duct flows and these include the electrical Hartmann number and electrical Reynolds number. The many studies communicated on either MHD duct or EHD duct flows have demonstrated that generally transverse magnetic field damps flows whereas axial (longitudinal) electrical field accelerates flows. These findings instigated interest in combining both electrical and magnetic fields in the science of electro-magnetohydrodynamics (EMHD) which emerged in the late 1960s. ¹⁹ The most comprehensive treatment of the mechanics of EMHD flows was presented by Eringen and Maugin ²⁰ who considered many types of fluids including ferromagnetic, dielectric, electro-viscous, ferroelectric, electro-strictive, electro-magneto-rheological (EMR), and hybrid magnetic ionic liquids. Subsequently many researchers have investigated the combined effects of magnetic and electrical fields (both static and alternating) in viscous duct flows. Jian and Chang ²¹ developed theoretical solutions for EMHD incompressible viscous flow in a horizontal a slit microchannel duct under a lateral uniform electrical field and a spatially non-uniform vertical magnetic field. They deployed the variation of parameter method (VPM), Gaussian quadrature and Chebyshev spectral collocation techniques, and obtained good correlation with experiments. They found that axial flow deceleration is induced with increasing magnetic field decay factor whereas strong acceleration is generated with larger electrical field strength parameter. Chakraborty and Paul ²² computed the EMHD transport in a microchannel, showing that at relatively weak magnetic field intensities, significant enhancement in volumetric flow rate is produced. They also noted that excessive electrical field can ramp up the Joule heating and induce wall degradation. Umavathi and Bég ²³ investigated the thermo-solutal EMHD nanofluid convection in a vertical duct for nuclear cooling applications. They showed that the polarity of the electric circuit strongly modifies the velocity distribution and interaction with magnetic field. For negative electrical field, the flow is directed downwards whereas for positive electrical field it is upwards. For the case of a short circuit (vanishing electrical field), a linear velocity distribution is computed. They also found that skin friction is depleted at the left duct wall with increasing Hartmann number and nanoparticle solid volume fraction but enhanced with those these parameters at the right wall. They also examined the relative performance of copper, titanium oxide, and silver nanoparticles and noted that silver nanoparticles generate the highest Nusselt number, whereas titanium oxide nanoparticles produce the lowest. Nusselt numbers were also shown to be suppressed with electrical field but enhanced with magnetic field.

The above studies ignored non-Newtonian effects. However, many studies have identified that in for example nuclear slurry duct flows and solar direct absorber tubes, the presence of particles suspended may lead to rheological behavior. This has further motivated researchers to explore the use of different non-Newtonian formulations to better characterize the nonlinear shear-stress-strain behavior of such fluids. Soundalgekar and Aranake ²⁴ used the Stokeian polar (couple stress) rheological model to study heat transfer in vertical duct flows. Muzychka and Edge ²⁵ developed an elegant model for power-law liquids based on the Rabinowitsch-Mooney formulation and computed volumetric flow rates in a variety of ducts including circular and elliptic tubes, parallel channels, rectangular ducts, isosceles triangular ducts, circular annular, and polygonal ducts. Other relevant studies include Chaudhuri et al. ²⁶ (on forced convection heat transfer in third-grade viscoelastic fluids in ducts) and Rajagopal and Na ²⁷ (third grade fluids with thermal buoyancy effects in ducts). All these studies confirmed the significant deviation in heat and momentum characteristics due to non-Newtonian behavior. A more sophisticated viscoelastic model of the differential type is the fourth-grade model. This is superior to the second grade and third grade Reiner-Rivlin models. In addition to predicting normal stress differences, the fourth-grade model can also accommodate the shear thinning or thickening characteristics of actual industrial slurries since it does not assume shear viscosity to be a constant). Several investigations have implemented this model in thermal duct flows. Kazemi et al. ²⁸ studied the hydromagnetic flow of a fourth-grade viscoelastic fluid in a duct with forced convection. They showed that there is a decrement in average velocity and the bulk temperature inside the duct with an increase in non-Newtonian parameter. Fatima et al. ²⁹ presented perturbation solutions for two-phase ionic fourth-grade flow in a convergent/divergent duct containing spherically homogeneous metallic particles under an axial electrical field. Nadeem and Ali ³⁰ conducted a second law thermodynamic analysis of fourth-grade viscoelastic flow with heat transfer in a vertical pipe using a homotopy method. Further investigations include Sajid et al. ³¹ , Salawu et al. ³² (on thermal stability of magnetized fourth-grade fluids), and Khan et al. (on bio-inspired pumping of chemically reacting fourth-grade fluids with species diffusion). These studies generally showed that with increasing viscoelastic material parameter associated with the fourth-grade model, deceleration is induced in the flow.

In nuclear and other energy duct applications, complex thermophysical phenomena arise. These include thermal dispersion, stratification, radiative heat transfer, and convective heating at the boundaries. Another interesting characteristic is quadratic convection (nonlinear convection or non-Boussinesq convection) which is associated with nonlinear density temperature differences in the buoyancy force term in natural convection. This generalizes the conventional Boussinesq approximation and provides more sophistication in thermal modeling of such systems. Temperature gradients can be significantly modified with quadratic effects and therefore considering the nonlinear variation of density with temperature is preferable to the simpler approach of a linear variation of density with temperature. A number of articles have addressed quadratic convection in multi-physical flows. Kunnegowda et al. ³³ used homotopy perturbation method to compute the MHD viscoplastic micro-duct transport with an exponential space dependent heat source and quadratic convection effects. Triveni and Mahanthesh ³⁴ simulated the nanofluid flow in an annular duct with quadratic convection and quadratic thermal radiation. Kumar and Sood ³⁵ studied the non-Boussinesq mixed convection stagnation point flow on a contracting surface to a porous medium permeated with a magnetic field. They showed that increasing quadratic convection parameter strongly modifies the transport and both accelerates the boundary layer and encourages thermal diffusion. Okoya et al. ³⁶ computed the natural convection flow from a moving wall with variable thermal conductivity and quadratic non-Boussinesq effects. Al-Kouz et al. ³⁷ evaluated quadratic convection and thermal radiation influences on second order viscoelastic boundary layer flow with a bvp5c finite difference algorithm in MATLAB. They showed that flow is accelerated with stronger quadratic convection and a cooling effect is induced. Further studies include Sabu et al. ³⁸ (on inclined plate MHD nanofluids), Bhatti et al. ³⁹ (on viscolastic Maxwell radiative-convection flow from a stretching wall with thermal relaxation effects), and Mahanthesh et al. ⁴⁰ (on dusty nanofluid convection from a vertical wall).

In the present work we theoretically study the fourth-grade EMHD natural convection flow in a vertical parallel plate duct with temperature-dependent viscosity and quadratic convection. Previous studies have considered alternative simpler rheological models including the Williamson model, ⁴¹ Phan-Thien-Tanner (PPT) model, ⁴² couple stress model, ⁴³ and second order Reiner-Rivlin model.^44,45 Furthermore, we incorporate the effects of viscous dissipation and Joule heating (Ohmic dissipation)^46,47 and the Reynolds exponential temperature-dependent viscosity model^48–52 all of which can exert marked modifications in nuclear duct electro-magneto-hydrodynamic flows. The novelty of the present work is therefore the collective consideration of fourth-grade viscoelasticity, viscous heating, exponential viscosity, and Joule dissipation, which have not been considered hitherto anywhere in the scientific literature, despite immediate applications in nuclear duct electromagnetic flow control. Fourth-grade fluids are materials with exceptional flow characteristics that fall between regular liquids and solids. These materials exhibit a complicated rheological behavior known as non-Newtonian flow, in which the viscosity (resistance to flow) varies with shear rate, temperature, and pressure. Fourth-grade fluid flow has applications in numerous scientific disciplines, including physics, chemistry, engineering, and material science. Understanding fourth-grade fluid is essential for comprehending the behavior of complex fluids such as emulsions, foams, and suspensions. Scientists can design sophisticated materials with the desired properties for a variety of applications, including drug delivery systems, cosmetics, and paint formulations, by studying the flow properties of these substances. Viscoelastic fourth-grade rheology is also applicable to geology and the study of geological formations. It explains phenomena such as lava flow, mudslides, and glacier movement. By comprehending non-Newtonian dynamics, scientists can predict and manage such natural occurrences more precisely. Fourth-grade fluid flow is crucial in a variety of other industries, including food processing, polymer manufacturing, and hydrocarbon extraction. Understanding the flow behavior of non-Newtonian substances is crucial for optimizing manufacturing processes, regulating product quality, and increasing productivity. In particular the modern thrust in electromagnetic smart polymer-based conducting materials in high-risk industries including nuclear reactor technology has mobilized great interest in electro-magneto-fluid dynamics of viscoelastic fourth-grade liquids. This provides a strong justification for the current study. The normalized nonlinear ordinary differential conservation equations and boundary conditions are solved with a homotopy perturbation technique. Detailed solutions for axial velocity, temperature, skin friction, and Nusselt number are presented graphically and in tables for the effects of key electromagnetic, thermal, and dissipation parameters. Detailed interpretation is included. The simulations as noted earlier are pertinent to providing a deeper understanding of laminar transport phenomena in nuclear and solar ducts deploying advanced functional rheological working liquids.

Mathematical model for EMHD non-Newtonian quadratic dissipative duct flow

The regime under investigation is visualized in Figure 1. An electromagnetic fourth-grade fluid flows under natural convection between the vertical parallel electrically insulated plates which are located a finite distance 2d apart. A Cartesian coordinate system ( $x\prime$ , $Y\prime$ ) is used. The walls of the duct are located at $x\prime =d$ and $x\prime =-d$ and are sustained at constant but non-similar temperatures ${\theta }_2$ and ${\theta }_1,$ where ${\theta }_1>{\theta }_2$ . The fourth-grade fluid is electrically conducting and incompressible. An axial electric field (E) and a transverse magnetic field ( $\stackrel{\rightharpoonup }{B}$ ) are imposed which are mutually orthogonal. Maxwell displacement, Hall current, and magnetic induction effects are ignored.

Figure 1.

Physical model for EMHD fourth-grade fluid flow in a vertical parallel duct.

Under these assumptions, the mass conservation (continuity) equation takes the form:

$$\nabla ·\overrightarrow{U}=0,$$

(1)

The momentum equation using Ohm’s law is defined as:

$$\begin{array}{c}\rho (\overrightarrow{U}·\nabla \overrightarrow{U}+\frac{\partial \overrightarrow{U}}{\partial t\prime })=-\nabla ·p\prime +\nabla \tau \prime +\stackrel{\rightharpoonup }{J}\times \stackrel{\rightharpoonup }{B}+g\hfill \\ [{\left(\rho \overline{\beta }\right)}_0(\theta -{\theta }_m)+{\left(\rho \overline{\beta }\right)}_1{(\theta -{\theta }_m)}^2],\hfill \end{array}$$

(2)

Here the stress tensor is designated by $\tau \prime$ , pressure is denoted by $p\prime ,$ g is gravity, linear and nonlinear thermal expansion coefficients are designated by ${\overline{\beta }}_0,{\overline{\beta }}_1$ , $\theta$ is the temperature, ${\theta }_m[=({\theta }_1+{\theta }_2)/2]$ is the mean temperature, density is designated by $\rho$ , time is $t\prime$ , and magnetic field is designated by $\stackrel{\rightharpoonup }{B}$ . The current density vector is designated by $\stackrel{\rightharpoonup }{J}$ and it is expressed as:

$$\overrightarrow{U}\times \stackrel{\rightharpoonup }{B}+\stackrel{\rightharpoonup }{E}),$$

(3)

Here electrical conductivity is designated by $\sigma$ , and electrical field strength is denoted by $\stackrel{\rightharpoonup }{E}$ .

The constitutive equation for fourth-grade Reiner-Rivlin differential viscoelastic fluid is expressed following ⁵³ as:

$$\begin{array}{c}\tau \prime =\sum _{i=1}^4{\mathbf{S}}_i,\\ {\mathbf{S}}_1=\mu \left(T\right){\hslash }_1,\\ {\mathbf{S}}_2={\alpha }_1{\hslash }_2+{\alpha }_2{\hslash }_1^2,\\ {\mathbf{S}}_3={\beta }_1{\hslash }_3+{\beta }_2\left({\hslash }_1{\hslash }_2+{\hslash }_2{\hslash }_1\right)+{\beta }_3\left(\mathrm{trace}{\hslash }_2\right){\hslash }_1,\\ \begin{array}{l}{\mathbf{S}}_4={\gamma }_1{\hslash }_4+{\gamma }_2\left({\hslash }_3{\hslash }_1+{\hslash }_1{\hslash }_3\right)+{\gamma }_3{\hslash }_2^2+{\gamma }_4\left({\hslash }_2{\hslash }_1^2+{\hslash }_1^2{\hslash }_2\right)\\ +{\gamma }_5\left(\mathrm{tr}{\hslash }_2\right){\hslash }_2+{\gamma }_6\left(\mathrm{tr}{\hslash }_2\right){\hslash }_1^2+\left[{\gamma }_7\mathrm{tr}{\hslash }_3+{\gamma }_8\mathrm{tr}{\hslash }_2{\hslash }_1\right]{\hslash }_1,\end{array}\end{array}\}$$

(4)

Here ${\alpha }_j,{\beta }_j,{\gamma }_j,j=1,2,3,...,8$ indicate the material constants, and the kinematical tensors ${\hslash }_1,{\hslash }_2,{\hslash }_3$ are expressed as:

$$\begin{array}{c}{\hslash }_1=\mathrm{\Xi }+{\mathrm{\Xi }}^t,\mathrm{\Xi }=\mathrm{grad}\overrightarrow{U},\hfill \\ {\hslash }_n=\frac{d{\hslash }_{n-1}}{dt\prime }+{\mathrm{\Xi }}^t{\hslash }_{n-1}+{\hslash }_{n-1}\mathrm{\Xi },n=2,3,\dots \hfill \end{array}\}$$

(5)

The energy equation considering the effects of Joule heating and viscous dissipation is specified as:

$$\left(\rho c_p\right)\frac{d\theta }{d{t}^{\prime }}=kgrad\theta +{\tau }^{\prime }:grad\overrightarrow{U}+\frac{\stackrel{\rightharpoonup }{J}·\stackrel{\rightharpoonup }{J}}{\sigma },$$

(6)

Here thermal conductivity of the fluid is designated by $k$ , $c_p$ designates the specific heat, and the symbols “·” and “:” represents the single and double dot vector products.

For the proposed problem, we seek the velocity field and the temperature field in the following form:

$$\overrightarrow{U}=[0,V(x\prime ),0],\theta =\theta (x\prime ).$$

(7)

The general vector equations may be reduced to the following component form:

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial p\prime }{\partial x\prime }=(2{\alpha }_1+{\alpha }_2)\frac{\mathrm{d}}{\mathrm{d}x\prime }{\left(\frac{\mathrm{d}V}{\mathrm{d}x\prime }\right)}^2}\\ +4({\gamma }_3+{\gamma }_4+{\gamma }_5+\frac{{\gamma }_6}{2})\frac{\mathrm{d}}{\mathrm{d}x\prime }{\left(\frac{\mathrm{d}V}{\mathrm{d}x\prime }\right)}^4,\hfill \end{array}$$

(8)

$$\begin{array}{c}\frac{\partial p″}{\partial y″}=\frac{\mathrm{d}}{\mathrm{d}x\prime }\left[\mu \left(\theta \right)\frac{\mathrm{d}V}{\mathrm{d}x\prime }\right]+6\left({\beta }_2+{\beta }_3\right)\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{x}^{\prime 2}}\left(\frac{\mathrm{d}V}{\mathrm{d}{x}^{\prime 2}}\right)-\sigma B^{2}V+\sigma EB\\ +\left[{\left(\rho \overline{\beta }\right)}_0\left(\theta -{\theta }_m\right)+{\left(\rho \overline{\beta }\right)}_1{\left(\theta -{\theta }_m\right)}^2\right]g,\end{array}\}$$

(9)

$$\frac{\partial p\prime }{\partial z\prime }=0$$

(10)

In equation (9), the viscoelastic fluid viscosity depends on temperature and the Reynolds exponential model is used which is defined as:

$$\mu \left(\theta \right)={\mu }_0\exp [-\overline{h}(\theta -{\theta }_m)],$$

(11)

Here ${\mu }_0$ indicates the dynamic viscosity, and $\overline{h}$ represents the strength dependency between $\theta$ and $\mu \left(\theta \right)$ .

The modified pressure is stated as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\overline{p}=p\prime -(2{\alpha }_1+{\alpha }_2){\left(\frac{\mathrm{d}V}{\mathrm{d}x\prime }\right)}^2}\\ -4({\gamma }_3+{\gamma }_4+{\gamma }_5+\frac{{\gamma }_6}{2}){\left(\frac{\mathrm{d}V}{\mathrm{d}x\prime }\right)}^4.\hfill \end{array}$$

(12)

In view of equation (12), equations (8–10) can be reduced to the following form as:

$$\frac{\partial \overline{p}}{\partial x\prime }=0,$$

(13)

$$\begin{array}{l}\frac{\partial \overline{p}}{\partial y\prime }=\frac{\mathrm{d}}{\mathrm{d}x″}\left[\mu \left(\theta \right)\frac{\mathrm{d}V}{\mathrm{d}x″}\right]+6\left({\beta }_2+{\beta }_3\right)\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{X}^{″2}}{\left(\frac{\mathrm{d}V}{\mathrm{d}X″}\right)}^2-\sigma B^{2}V+\sigma EB\\ +\left[{\left(\rho \overline{\beta }\right)}_0\left(\theta -{\theta }_m\right)+{\left(\rho \overline{\beta }\right)}_1{\left(\theta -{\theta }_m\right)}^2\right]g,\end{array}\}$$

(14)

$$\frac{\partial \overline{p}}{\partial z\prime }=0.$$

(15)

The preceding equation indicates that the pressure gradient $\partial \hat{p}/\partial y\prime$ is constant. We obtain:

$$\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}x″}\left[\mu \left(\theta \right)\frac{\mathrm{d}V}{\mathrm{d}x″}\right]+6\left({\beta }_2+{\beta }_3\right)\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{x}^{\prime 2}}{\left(\frac{\mathrm{d}V}{\mathrm{d}x\prime }\right)}^2-\sigma B^{2}V+\sigma EB\\ +\left[{\left(\rho \beta \right)}_0\left(\theta -{\theta }_m\right)+{\left(\rho \beta \right)}_1{\left(\theta -{\theta }_m\right)}^2\right]g=\frac{\partial \overline{p}}{\partial y\prime },\end{array}\}$$

(16)

The following are the velocity boundary conditions at the duct walls:

$$V(\pm d)=0,$$

(17)

The energy equation including viscous dissipation and Joule heating is defined as:

$$\begin{array}{c}k\frac{{\mathrm{d}}^2\theta }{\mathrm{d}x{\prime }^2}+\mu \left(\theta \right){\left(\frac{\mathrm{d}V}{\mathrm{d}x{\prime }^2}\right)}^2+2({\beta }_2+{\beta }_3){\left(\frac{\mathrm{d}V}{\mathrm{d}x\prime }\right)}^4\hfill \\ +\sigma (B^2{V}^2-2\mathrm{BEV}+E^2)=0.\hfill \end{array}$$

(18)

The relevant thermal boundary conditions prescribed at the duct walls are as follows:

$$\theta (-d)={\theta }_1,\theta (+d)={\theta }_2,$$

(19)

The next step is to transform the above equations into dimensionless form using the scaling variables listed below:

$$v=\frac{V}{v_0},x=\frac{x\prime }{d},T=\frac{\theta -{\theta }_m}{{\theta }_1-{\theta }_2},$$

(20)

Here $v_0$ indicates the reference velocity. Using equation (20), the dimensionless form of equations (16) and (19) emerge as:

$$\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}x}\left(\exp [-\epsilon T]\frac{\mathrm{d}v}{\mathrm{d}x}\right)+6\chi {\left(\frac{\mathrm{d}v}{\mathrm{d}x}\right)}^2\frac{{\mathrm{d}}^{2}v}{\mathrm{d}{x}^2}-{\gamma }^{2}v\hfill \\ +\eta +{\lambda }_r(T+{\lambda }_c{T}^2)={\mathrm{\Omega }}_p,\hfill \end{array}$$

(21)

$$\frac{{\mathrm{d}}^{2}T}{\mathrm{d}{x}^2}+\delta \exp [-\epsilon T]{\left(\frac{\mathrm{d}v}{\mathrm{d}x}\right)}^2+2\delta \chi {\left(\frac{\mathrm{d}v}{\mathrm{d}x}\right)}^4+\delta {\gamma }^2{v}^2-{\zeta }_{1}v+{\zeta }_2=0.$$

(22)

The boundary conditions (17) and (19) now assume the dimensionless form:

$$\begin{array}{c}v(-1)=0,v(+1)=0,\hfill \\ T(-1)=0.5,T(+1)=-0.5,\hfill \end{array}\}$$

(23)

And

$$\begin{array}{c}\chi =\frac{v_0^2\left({\beta }_2+{\beta }_3\right)}{{\mu }_0{d}^2},\mathrm{Re}=\frac{v_{0}d}{\upsilon },Gr=\frac{g{d}^3{\left(\rho \overline{\beta }\right)}_0\left({\theta }_1-{\theta }_2\right)}{\upsilon {\mu }_0},\\ \epsilon =\overline{h}\left({\theta }_1-{\theta }_2\right),\eta =\frac{\sigma BE{d}^2}{{\mu }_0{v}_0},{\lambda }_c=\frac{{\left(\rho \overline{\beta }\right)}_0\left({\theta }_1-{\theta }_2\right)}{{\left(\rho \beta \right)}_1},\\ \gamma =Bd\sqrt{\frac{\sigma }{{\mu }_0}},{\lambda }_r=\frac{Gr}{\mathrm{Re}},{\Omega }_p=\frac{d^2}{v_0{\mu }_0}\frac{\partial \overline{p}}{\partial y\prime }\delta =\frac{{\mu }_0{v}_0^2}{k\left({\theta }_1-{\theta }_2\right)},\\ {\zeta }_1=\frac{2BE{d}^2{v}_0\sigma }{\left({\theta }_1-{\theta }_2\right)k},{\zeta }_2=\frac{d^2{E}^2\sigma }{\left({\theta }_1-{\theta }_2\right)k},\upsilon =\frac{{\mu }_0}{\rho }.\end{array}\}$$

(24)

Here $\gamma$ is Hartmann number, Re is Reynolds number, Gr is the thermal Grashof number, $\epsilon$ indicates the variable viscosity parameter, ${\mathrm{\Omega }}_p$ is the dimensionless pressure gradient, $\eta$ is electrical field strength parameter, ${\lambda }_c$ quadratic thermal convection parameter, $\chi$ is fourth-grade fluid parameter, ${\lambda }_r$ is thermal buoyancy parameter, $\delta$ is Brinkman number (which indicates the ratio of heat created by viscous dissipation to heat transferred by molecular conduction), ${\zeta }_1$ is a parameter relating the effects of heat generation owing to the interaction of magnetic and electric fields on heat conduction, ${\zeta }_2$ is ratio of Joule heating to heat conduction, and ν is kinematic viscosity of the fourth-grade fluid.

On the left plate $x\prime =-d$ , the skin friction coefficient and Nusselt number are expressed in dimensionless format following Rajagopal and Na ²⁷ :

$$\begin{array}{c}S~{\exp [-\epsilon T]\frac{\mathrm{dv}}{\mathrm{dx}}|}_{x=-1}+{2\gamma {\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^3|}_{x=-1},\hfill \\ h~{-\frac{\mathrm{dT}}{\mathrm{dx}}|}_{x=-1}.\hfill \end{array}\}$$

(25)

Perturbation solutions using HPM

The transformed ordinary differential equations for momentum and energy, that is, equations (21) and (22) are nonlinear and coupled, rendering exact solutions intractable. As a result, we utilize the homotopy perturbation method (HPM) for determining series solutions. This method was introduced as a special case of the homotopy analysis method (HAM) and is very adaptive and yields excellent accuracy and stability for nonlinear differential equation systems. It has been utilized in a wide spectrum of complex engineering science problems, including structural vibration, ⁵⁴ peristaltic pumping,^55,56 and more recently, electromagnetic duct flows containing porous media. ⁴⁴ The perturbation formulation for the present problem is as follows:

$$\begin{array}{l}h\left(u,\zeta \right)=\left(1-\zeta \right)\left[L_0\left(u\right)-L_0\left({\overline{v}}_0\right)\right]+\zeta \\ [L_0\left(u\right)-\epsilon \mathrm{\Theta }\frac{{\mathrm{d}}^{2}u}{\mathrm{d}{x}^2}-\epsilon \frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}x}\frac{\mathrm{d}u}{\mathrm{d}x}+6\chi {\left(\frac{\mathrm{d}u}{\mathrm{d}x}\right)}^2\frac{{\mathrm{d}}^{2}u}{\mathrm{d}{x}^2}-{\gamma }^{2}u+\eta +{\lambda }_r\left(\mathrm{\Theta }+{\lambda }_c{\mathrm{\Theta }}^2\right)-{\Omega }_p],\end{array}$$

(26)

$$\begin{array}{l}h\left(\mathrm{\Theta },\zeta \right)=\left(1-\zeta \right)\left[L_0\left(\mathrm{\Theta }\right)-L_0\left({\overline{T}}_0\right)\right]+\zeta \\ [L_0\left(\mathrm{\Theta }\right)+\delta \left(1-\epsilon \mathrm{\Theta }\right){\left(\frac{\mathrm{d}u}{\mathrm{d}x}\right)}^2+2\delta \chi {\left(\frac{\mathrm{d}u}{\mathrm{d}x}\right)}^4+\left(\delta {\gamma }^2{u}^2-{\zeta }_{1}u+{\zeta }_2\right)],\end{array}$$

(27)

Where $\zeta$ is the embedding parameter. The linear operator $L_0$ is selected and the following tentative guesses, $u_0,{\mathrm{\Theta }}_0$ are made:

$$\begin{array}{c}{L}_0=\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}\},\mathrm{Linear}\mathrm{operator}\hfill \\ {\overline{v}}_0=\frac{1}{4}(x^2-1),{\mathrm{\Theta }}_0=-\frac{x}{2}\},\mathrm{Initial}\mathrm{guess}\hfill \end{array}$$

(28)

Let us now propose the series expansions for equations (26) and (27), that is,

$$\{\begin{array}{c}\multicolumn{1}{c}{u=u_0+\zeta u_1+{\zeta }^2{u}_2+\dots ,}\\ \multicolumn{1}{c}{\theta ={\theta }_0+\zeta {\theta }_1+{\zeta }^2{\theta }_2+\dots .}\end{array}$$

(29)

By inserting equation (29) into equations (26) and (27), we acquire a set of linear differential equations at every level. After the accomplishment of this stage, we will utilize the symbolic software Mathematica to determine the linear differential equation solutions for each order.

Zeroth other system ${\zeta }^0$ and solutions

By substituting equation (29) into equations (26) and (27), we get the zeroth order system as follows:

$$\begin{array}{c}{L}_0\left(u_0\right)-L_0\left({\overline{v}}_0\right)=0,\hfill \\ u_0(\pm 1)=0,\hfill \end{array}\}$$

(30)

$$\begin{array}{c}{L}_0\left({\mathrm{\Theta }}_0\right)-L_0\left({\overline{T}}_0\right)=0,\hfill \\ {\overline{T}}_0(-1)=0.5,{\overline{T}}_0(+1)=-0.5,\hfill \end{array}\}$$

(31)

The zeroth order solution is found to be:

$$\begin{array}{c}{u}_0=\frac{x^2-1}{4},\hfill \\ {\mathrm{\Theta }}_0=-\frac{x}{2},\hfill \end{array}\}$$

(32)

First order system ${\zeta }^1$ and solutions

The first order system is obtained in the following form:

$$\begin{array}{c}{L}_0\left(u_1\right)+L_0\left({\overline{v}}_0\right)-\epsilon {\mathrm{\Theta }}_0\frac{{\mathrm{d}}^2{u}_0}{\mathrm{d}{x}^2}-\epsilon \frac{\mathrm{d}{\mathrm{\Theta }}_0}{\mathrm{d}x}\frac{\mathrm{d}{u}_0}{\mathrm{d}x}+6\chi {\left(\frac{\mathrm{d}{u}_0}{\mathrm{d}x}\right)}^2\frac{{\mathrm{d}}^2{u}_0}{\mathrm{d}{x}^2}\hfill \\ -{\gamma }^2{u}_0+\eta +{\lambda }_r({\mathrm{\Theta }}_0+{\lambda }_c{\mathrm{\Theta }}_0^2)-{\mathrm{\Omega }}_p,\hfill \\ u_1(\pm 1)=0,\hfill \end{array}\}$$

(33)

$$\begin{array}{c}{L}_0\left({\mathrm{\Theta }}_1\right)+L_0\left({\overline{T}}_0\right)+\delta (1-\epsilon {\mathrm{\Theta }}_0){\left(\frac{\mathrm{d}{u}_0}{\mathrm{d}x}\right)}^2+2\delta \chi {\left(\frac{\mathrm{d}{u}_0}{\mathrm{d}x}\right)}^4\hfill \\ +(\delta {\gamma }^2{u_0}^2-{\zeta }_1{u}_0+{\zeta }_2),\hfill \\ {\mathrm{\Theta }}_1(\pm 1)=0,\hfill \end{array}\}$$

(34)

The solution of the first order system is obtained as:

$$\begin{array}{c}\multicolumn{1}{c}{u_1=\frac{x^2-1}{48}[H{a}^2(-5+x^2)-{\lambda }_c{\lambda }_r+4x(-\epsilon +\lambda r)}\\ -x^2({\lambda }_c{\lambda }_r+3\chi )-3(4+8\eta +\chi -8{\mathrm{\Omega }}_p)],\hfill \end{array}$$

(35)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{\Theta }}_1=\frac{1-x^2}{480}[{\gamma }^2(11-4x^2+x^4)\delta -10(-5+x^2){\zeta }_1}\\ +240{\zeta }_2+2\delta \{(1+x^2)(5+\chi )+x^4\chi \}],\hfill \end{array}$$

(36)

Second order system ${\zeta }^2$ and solutions

The second order system is obtained as:

$$\begin{array}{c}{L}_0\left(u_2\right)-\epsilon {\mathrm{\Theta }}_1\frac{{\mathrm{d}}^2{u}_0}{\mathrm{d}{x}^2}-\epsilon {\mathrm{\Theta }}_0\frac{{\mathrm{d}}^2{u}_1}{\mathrm{d}{x}^2}-\epsilon \frac{\mathrm{d}{\mathrm{\Theta }}_1}{\mathrm{d}x}\frac{\mathrm{d}{u}_0}{\mathrm{d}x}-\epsilon \frac{\mathrm{d}{\mathrm{\Theta }}_0}{\mathrm{d}x}\frac{\mathrm{d}{u}_1}{\mathrm{d}x}\hfill \\ +6\chi [2\frac{\mathrm{d}{\overline{v}}_0}{\mathrm{d}x}\frac{{\mathrm{d}}^2{\overline{v}}_0}{\mathrm{d}{x}^2}\frac{\mathrm{d}{\overline{v}}_1}{\mathrm{d}x}+\frac{{\mathrm{d}}^2{\overline{v}}_1}{\mathrm{d}{x}^2}{\left(\frac{\mathrm{d}{\overline{v}}_0}{\mathrm{d}x}\right)}^2]\hfill \\ -{\gamma }^2{u}_1+{\lambda }_r({\mathrm{\Theta }}_1+2{\lambda }_c{\mathrm{\Theta }}_0{\mathrm{\Theta }}_1),\hfill \\ u_2(\pm 1)=0,\hfill \end{array}\}$$

(37)

$$\begin{array}{c}{L}_0\left({\mathrm{\Theta }}_1\right)+\delta [2\frac{\mathrm{d}{u}_1}{\mathrm{d}x}+8\chi \frac{\mathrm{d}{u}_1}{\mathrm{d}x}{\left(\frac{\mathrm{d}{u}_0}{\mathrm{d}x}\right)}^2-\epsilon {\mathrm{\Theta }}_1-2\epsilon {\mathrm{\Theta }}_0\frac{\mathrm{d}{u}_1}{\mathrm{d}x}]\hfill \\ \frac{\mathrm{d}{u}_0}{\mathrm{d}x}+2\delta {\gamma }^2{u}_1{u}_0-{\zeta }_1{u}_1,\hfill \\ {\mathrm{\Theta }}_1(\pm 1)=0,\hfill \end{array}\}$$

(38)

The solution of the second order system is obtained as:

$$\begin{array}{c}\multicolumn{1}{c}{u_2=u_{2,1}+u_{2,2}{x}^1+u_{2,3}{x}^2+u_{2,4}{x}^3+u_{2,5}{x}^4+u_{2,6}{x}^5}\\ +u_{2,7}{x}^6+u_{2,8}{x}^7+u_{2,9}{x}^8+u_{2,10}{x}^9,\hfill \end{array}$$

(39)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{\Theta }}_2={\mathrm{\Theta }}_{2,1}+{\mathrm{\Theta }}_{2,2}{x}^1+{\mathrm{\Theta }}_{2,3}{x}^2+{\mathrm{\Theta }}_{2,4}{x}^3+{\mathrm{\Theta }}_{2,5}{x}^4}\\ +{\mathrm{\Theta }}_{2,6}{x}^5+{\mathrm{\Theta }}_{2,7}{x}^6+{\mathrm{\Theta }}_{2,8}{x}^7+{\mathrm{\Theta }}_{2,9}{x}^8+{\mathrm{\Theta }}_{2,10}{x}^9.\hfill \end{array}$$

(40)

The constants $u_{2,n},{\mathrm{\Theta }}_{2,n};n=1,2,\dots 10,$ derived in the preceding calculations can be found in the Appendix. We may deduce the final form of the outcomes employing the HPM condition as:

$$v=\underset{\zeta \to 1}{lim}u=u_0+\zeta u_1+{\zeta }^2{u}_2+\dots .$$

(41)

$$\theta =\underset{\zeta \to 1}{lim}{\theta }_0+\zeta {\theta }_1+{\zeta }^2{\theta }_2+\dots .$$

(42)

Results and discussion

In this section, we present detailed graphical and tabular solutions for velocity, temperature, skin friction, and the Nusselt number in the duct, for the influence of all key control parameters. All computations are executed in Mathematica and shown in Table 1 and Figures 2 to 15. As one particular case of our investigation, Table 2 displays the numerical comparison with previously published data. It is noticed that the current findings are precisely consistent with Abbasi et al. ⁵⁷ .

Table 1.

Comparison of the velocity profile with previously published results by fixing the remaining parameters as $\epsilon =\gamma =\eta ={\lambda }_r={\lambda }_c=0,\chi =0.5,{\mathrm{\Omega }}_p=1$ .

x	Present results	Abbasi et al. 57
0	−0.395353044732200	−0.395353045
0.1	−0.390278044901173	−0.390377559
0.2	−0.375724400825540	−0.375724401
0.3	−0.352088141541283	−0.352088141
0.4	−0.320513599352345	−0.320336095
0.5	−0.281320880996496	−0.281320881
0.6	−0.235797999318292	−0.235797999
0.7	−0.184408971096291	−0.184408971
0.8	−0.127691231510506	−0.127691231
0.9	−0.066095053134708	−6.61E-02
1	0	0

Figure 2.

Coefficient of skin friction for multiple values of variable viscosity parameter $\epsilon$ , fourth-grade fluid parameter $\chi$ , Brinkman number $\delta$ , electric field strength parameter $\eta$ , Hartmann number $\gamma$ , thermal buoyancy parameter ${\lambda }_r$ , and quadratic thermal convection parameter ${\lambda }_c$ .

Figure 3.

Nusselt number for multiple values of variable viscosity parameter $\epsilon$ , fourth-grade fluid parameter $\chi$ , Brinkman number $\delta$ , electric field strength parameter $\eta$ , Hartmann number $\gamma$ , thermal buoyancy parameter ${\lambda }_r$ , and quadratic thermal convection parameter ${\lambda }_c$ .

Figure 4.

Influence of fourth-grade material parameter $\chi$ on velocity distribution.

Figure 5.

Influence of electric field strength parameter $\eta$ on velocity distribution.

Figure 6.

Influence of variable viscosity parameter $\epsilon$ on velocity distribution.

Figure 7.

Influence of thermal buoyancy (natural) convection parameter ${\lambda }_r$ on velocity distribution.

Figure 8.

Effect of quadratic thermal convection parameter ${\lambda }_c$ on velocity distribution.

Figure 9.

Influence of Hartmann (magnetic) number $\gamma$ on velocity distribution.

Figure 10.

Effect of fourth-grade material fluid parameter $\chi$ on temperature distribution.

Figure 11.

Influence of Brinkman number $\delta$ on temperature distribution.

Figure 12.

Influence of electromagnetic heat generation to conduction ratio ${\zeta }_1$ on temperature distribution.

Figure 13.

Influence of variable viscosity parameter $\epsilon$ on temperature distribution.

Figure 14.

Effect of thermal buoyancy parameter ${\lambda }_r$ on temperature distribution.

Figure 15.

Effect of quadratic thermal convection parameter ${\lambda }_c$ on temperature distribution.

Table 2.

HPM results for the effects of key parameters on skin friction and Nusselt number.

Parameters	Values	Coefficient of skin friction	Nusselt number
$\mathit{\epsilon }$	0	0.363579642	0.058035341
	0.4	0.346855352	0.064830629
	0.6	0.327962242	0.073336828
$\mathit{\chi }$	0	0.335917849	0.070981940
	0.5	0.343657214	0.066632321
	1	0.348663956	0.063809635
$\mathit{\delta }$	0.3	0.345738770	0.070702503
	0.5	0.346855352	0.064830629
	0.7	0.347993394	0.058866085
$\mathit{\eta }$	1	0.346855352	0.064830629
	1.5	0.776411238	0.110749135
	2	1.240973934	0.096432223
$\mathit{\gamma }$	0	0.363702482	0.069451489
	0.5	0.346855352	0.064830629
	1	0.306876751	0.053988768
${\mathit{\lambda }}_{\mathit{r}}$	0.5	0.176197991	0.040934872
	1	0.346855352	0.064830629
	2	0.685748678	0.082449520
${\mathit{\lambda }}_{\mathit{c}}$	0	0.285794633	0.057846555
	0.5	0.346855352	0.064830629
	1	0.407345992	0.070554752
${\mathit{\zeta }}_1$	0	0.362327892	0.025694391
	0.5	0.346855352	0.064830629
	1	0.333674090	0.099221419

Skin friction and Nusselt number

Table 2 shows skin friction coefficient and Nusselt number for all the emerging parameters.

A significant enhancement in skin friction accompanies an elevation in fourth-grade fluid parameter $\chi$ , Brinkman number $\delta$ , electric field parameter $\eta$ , thermal buoyancy parameter ${\lambda }_r,$ and quadratic thermal convection parameter ${\lambda }_c$ . However, skin friction is strongly reduced with a rise in variable viscosity parameter $\epsilon$ , Hartmann (magnetic) number $\gamma$ , and electromagnetic heat generation to conduction ratio, ${\zeta }_1$ . Nusselt number magnitudes are elevated with increase in variable viscosity parameter $\epsilon$ , electromagnetic heat generation to conduction ratio, ${\zeta }_1$ thermal buoyancy parameter ${\lambda }_r$ , and quadratic thermal convection parameter ${\lambda }_c$ , whereas they are significantly suppressed with increment in fourth-grade fluid parameter $\chi$ , Brinkman number $\delta$ , and Hartmann magnetic number $\gamma .$ A more complex response is computed for increment in electric field parameter $\eta$ − initially the Nusselt number is boosted from $\eta$  = 0 (vanishing electrical field) to $\eta$  = 1, but a subsequent increase to $\eta$  = 1.5 actually decreases heat transfer to the wall of the duct.

The results In Table 2 are further visualized in Figures 2 and 3. Figure 2 shows the variation in the skin friction profile. It is noticed here that increment in the variable viscosity parameter $\epsilon$ and the Hartmann magnetic field produce a consistent decrement in skin friction profile. Clearly a boost in exponential viscosity will enhance overall viscosity of the fourth-grade fluid and this will inhibit shearing along the wall, manifesting in a lower skin friction. A stronger magnetic field will also oppose momentum transfer and will decelerate the flow reducing the skin friction at the wall but simultaneously enhancing the momentum boundary layer thickness. The fourth-grade material parameter has a different influence on skin friction since it is related to stress relaxation and retardation, not merely viscosity. An increase in this parameter causes a gradual increase in skin friction implying a thinner momentum boundary layer. An increment followed by a weak increase is observed in skin friction with Brinkman number due to kinetic energy dissipation (internal friction). Increasing thermal buoyancy parameter strongly elevates the skin friction due to the intensity of natural convection currents which assist momentum development. This assists the flow and manifests in an elevation in velocity which implies faster shearing of the fluid against the duct walls and produces the associated hike in skin friction. Similarly increasing electrical field strength parameter and quadratic convection both enhance skin friction and produce thinner momentum boundary layers. Clearly the stronger electrical field intensity assists axial momentum development and is the primary reason for deploying such an orientation of the electrical field. The nonlinear effect of quadratic convection therefore avoids the under-prediction in skin friction associated with the classical Boussinesq (linear) model.

Figure 3 shows that Nusselt number is elevated with variable viscosity parameter but is suppressed with fourth-grade material parameter. A significant decrement in Nusselt number is also associated with greater Brinkman number since temperatures are elevated in the boundary layer and the net effect is to draw heat away from the wall, leading to a cooling there. A thinner thermal boundary layer thickness will result with higher Brinkman numbers associated with stronger viscous dissipation. Hartmann number elevation also induces a strong downturn in Nusselt number since the electromagnetic fluid is heated with higher magnetic field. Increasing the magnetic field in the fluid flow amplifies Lorentz force that dominates over the viscous force. This results in a net transfer of heat away from the wall into the core region of the duct and decreases Nusselt number. A similar though less prominent increment in Nusselt number accompanies an increase in thermal buoyancy and quadratic convection parameters. Again, the initial increase in Nusselt number is computed with a rise in electrical field strength parameter but subsequent increase in this parameter produces the opposite trend which is associated with the presence of Joule dissipation in the flow.

Velocity distribution

Figures 4 to 9 show the impact of all key parameters on velocity profiles. The velocity profile in Figure 4 exhibits a strong depletion with larger values of the fourth-grade fluid parameter. In practical terms, an increase in $\chi$ increases the stress relaxation relative to retardation in a fourth-grade which counteracts momentum diffusion and decelerates the flow. This will increase hydrodynamic boundary layer thickness at the duct walls. Furthermore, when $\chi =0$ , the classical case of a Newtonian viscous fluid is retrieved. Figure 5 indicates the velocity profile is significantly enhanced with greater electric field intensity $(\eta )$ . When an electric field is introduced into a working fluid containing charged particles, the particles are attracted to the electrodes, causing the fluid to be induced in the same direction as the applied electric field. This technique is most often seen in microfluidic devices when a strong electric field is present in relatively small channels. However, it can also be used in larger nuclear reactor ducts to regulate momentum distribution. Strong acceleration is therefore produced with axial electrical field. Figure 6 depicts the effect of the varying viscosity parameter $\epsilon$ on the velocity profile. It is evident that the velocity increases in the left half space but then decreases in the right half space with increment in viscosity parameter. The peak velocity arises consistently near the centerline of the duct. Asymmetrical parabolic profiles span the duct width. Figure 7 depicts the effect of the thermal buoyancy convection parameter ${\lambda }_r$ on the velocity profile. With increasing values of this parameter, the flow is accelerated across the channel width, that is, velocity profile is greatly enhanced owing to the enhanced contribution of thermal buoyancy relative to viscous forces. For ${\lambda }_r$  = 0 forced convection is present, and velocities are minimized. For ${\lambda }_r$  = 1 both thermal buoyancy and viscous force contribute equally. When ${\lambda }_r$  > 1, thermal buoyancy dominates the regime and significant momentum flux escalation is produced. It is also noteworthy that as viscosity varies with temperature and shear rate, the fluid becomes less dense in the presence of heat and this induces a boost to momentum resulting in a hike in velocity. Figure 8 demonstrates that acceleration is also produced in the axial flow with increasing quadratic convection parameter ${\lambda }_c$ . The nonlinear temperature dependence of the non-Boussinesq buoyancy term directly contributes to the elevation in velocity. Figure 9 depicts the impact of the Hartmann number $\gamma$ on the velocity field. As can be observed, the velocity profile decreases consistently owing to the significant retarding influence of the magnetic field. When a magnetic field is applied to a working fluid, it experiences a Lorentz force, which is a drag force that opposes fluid motion. Strong flow deceleration is therefore induced with higher Hartmann number across the entire span of the duct. Clearly the non-magnetic scenario ( $\gamma$  = 0), that is, vanishing magnetic field produces the maximum velocity. For $\gamma$  = both Lorentz magnetic drag force and viscous hydrodynamic force contribute equally. For $\gamma$  > 1 the Lorentz magnetic drag force dominates the viscous force leading to strong flow deceleration. Symmetric parabolic topologies are computed across the duct. The excellent flow regulation properties of transverse magnetic field are thereby confirmed.

Temperature distribution

Figures 10 to 15 depict the evolution of the temperature profiles with key parameters. Figure 10 demonstrates that the temperature is boosted with elevation in fourth-grade fluid $\chi$ parameter rises. However, we can see that the implications are insignificant closer to the channel’s walls, with only a minor increase in the middle of the channel. The heating of the fluid is associated with stress relaxation and retardation effects. Thermal diffusion is encouraged with these phenomena. Figure 11 illustrates that enhancing the Brinkman number $\delta$ causes the temperature profile to strongly increase, in particular, near the left wall as dissipation enhances the conversion of kinetic energy to heat via internal friction. As a result, the fourth-grade fluid also decelerates. The case of vanishing Brinkman number, $\delta$  = 0 clearly corresponds to the minimal temperature and confirms that non-inclusion of viscous heating leads to an under-prediction in temperature in the bulk fluid. For this case the fourth order derivative term, $2\delta \gamma {\left(\frac{\mathrm{dv}}{\mathrm{dx}}\right)}^4$ in equation (22) disappears. Figure 12 demonstrates the fact that the electromagnetic heat generation relative to thermal conduction parameter ${\zeta }_1$ strongly suppresses temperature profile. Clearly when this parameter vanishes ${\zeta }_1$  → 0, and thermal conduction dominates over the electrical and magnetic field heating contribution and maximum temperatures are computed. Figure 13 shows that the variable viscosity parameter produces a opposite behavior in temperature, closer to the duct centerline. The exponential viscosity variation exerts a much less dramatic effect on temperature than the fourth-grade material parameter. Figure 14 shows that there is a slight decrement in temperatures across the width of the duct with increment in thermal buoyancy (natural convection) parameter, ${\lambda }_r$ . As with the other temperature profiles, the weak temperature is always computed at the left hotter duct wall and decays monotonically to the right cooler wall. Thermal buoyancy therefore primarily influences the momentum field, and indeed features only in that equation. Via coupling with the energy (heat) equation there is an indirect effect on the temperature distribution. The temperature, however, is not as markedly influenced as the velocity field. Similarly in Figure 15 there is a relatively weak decrease in temperature across the duct width with greater values of quadratic convection parameter ${\lambda }_r$ . Again, the minor decrement is associated with the coupling effect between the momentum and energy equations (21) and (22), specifically via the term, ${\lambda }_r(T+{\lambda }_c{T}^2)$ in equation (21). Generally smooth profiles are computed for both temperature and velocity distributions testifying to the excellent numerical stability of HPM.

Concluding remarks

A new theoretical study has been presented for fourth-grade viscoelastic EMHD natural convection flow in a vertical parallel plate duct with non-Boussinesq quadratic convection. The effects of viscous dissipation and Joule heating (Ohmic dissipation) and Reynolds exponential temperature-dependent viscosity model have been included, all of which arise in nuclear and solar hybrid duct transport phenomena. The normalized nonlinear ordinary differential boundary value problem has been solved with a homotopy perturbation method (HPM) up to second order approximation which achieves exceptional accuracy and stability in terms of power series solutions. Detailed solutions for axial velocity, temperature, skin friction, and Nusselt number have been depicted graphically and in tables for the effects of key electromagnetic, thermal, and dissipation parameters. The key findings of the current analysis may be summarized as follows:

• (i) There is a considerable enhancement in wall skin friction with an increment in fourth-grade fluid parameter $\chi$ , Brinkman number $\delta$ , electrical field parameter $\eta$ , thermal buoyancy parameter ${\lambda }_r$ , and quadratic thermal convection parameter ${\lambda }_c$ . However, skin friction is strongly reduced with a rise in variable viscosity parameter $\epsilon$ , Hartmann (magnetic) number $\gamma$ , and electromagnetic heat generation to conduction ratio, ${\zeta }_1$ .

• (ii) Nusselt number magnitudes are elevated with increase in variable viscosity parameter $\epsilon$ and electromagnetic heat generation to conduction ratio, ${\zeta }_1$ , thermal buoyancy parameter ${\lambda }_r$ , and quadratic thermal convection parameter ${\lambda }_c,$ whereas they are significantly suppressed with increment in fourth-grade fluid parameter $\chi$ , Brinkman number $\delta$ , and Hartmann magnetic number $\gamma$ . The latter parameters therefore achieve better thermal control at the duct walls which is critical in regulating high temperatures (“thermal management”) in real nuclear reactor duct transport.

• (iii) A variable response in Nusselt number is computed for increment in electric field parameter $\eta ;$ initially the Nusselt number is boosted from $\eta$  = 0 (vanishing electrical field) to $\eta$  = 1, but a subsequent increase to $\eta$  = 1.5 actually decreases heat transfer to the wall of the duct.

• (iv) Increasing fourth-grade viscoelastic parameter $\chi$ elevates the viscosity and damps the axial velocity strongly as does higher magnetic Hartmann number.

• (v) A minor increase in temperature is induced with greater viscoelastic parameter $\chi$ , while a minor decrement is noticed for thermal buoyancy parameter ${\lambda }_r$ and quadratic thermal convection parameter ${\lambda }_c$ .

• (vi) A marked boost in temperature is generated with increasing Brinkman number $\delta$ whereas a very minor reduction is produced with a large increment in variable viscosity parameter.

• (vii) Increasing electrical strength parameter $\eta$ produces significant flow acceleration.

• (viii) In the left half space of the duct an increase in variable viscosity parameter $\epsilon$ produces significant axial flow acceleration whereas in the right half space weak deceleration is produced.

• (ix) Consistent acceleration, that is, higher velocities are associated with an increment in thermal buoyancy parameter ${\lambda }_r$ (ratio of thermal Grashof number and Reynolds number).

• (x) A weaker but distinct increase in axial velocity is also observed with higher values of quadratic thermal convection parameter ${\lambda }_c$ .

• (xi) The greatest elevation in skin friction number is produced with high electrical field strength parameter whereas the largest Nusselt number at the wall corresponds to low Brinkman number.

• The current study has probed deeper into emerging smart electromagnetic liquids for nuclear and other energy systems duct transport phenomena. Attention has, however, been confined to steady state flows. Future work may consider time-dependent effects and additionally the use of other non-Newtonian microstructural models such as the polar and micropolar family of formulations. Furthermore, alternating electrical and magnetic fields may also be addressed to reflect further industrial scenarios.