Irreversibility Analysis for Eyring–Powell Nanoliquid Flow Past Magnetized Riga Device with Nonlinear Thermal Radiation

Abstract

The report contained in this article is based on entropy generation for a reactive Eyring–Powell nanoliquid transfer past a porous vertical Riga device. In the developed model, the impacts of viscous dissipation, thermophoresis alongside nonlinear heat radiation and varying heat conductivity are modelled into the heat equation. The dimensionless transport equations are analytically tackled via Homotopy analysis method while the computational values of chosen parameters are compared with the Galerkin weighted residual method. Graphical information of the various parameters that emerged from the model are obtained and deliberated effectively. The consequences of this study are that the temperature field expands with thermophoresis, Brownian motion and temperature ratio parameters as the modified Hartmann number compels a rise in the velocity profile. The entropy generation rises with an uplift in fluid material term as well as Biot and Eckert numbers whereas Bejan number lessens with Darcy and Eckert parameters.

1. Introduction

Recently, a substantial amount of time and energy have been dedicated to examine non-Newtonian liquids because of their many uses in food manufacturing, polymers, chemical processing industry, crude oil extrusion, and so on. With different types of geometries and constitutive models, the mathematics of non-Newtonian fluids are non-linearly formulated, which are often difficult to solve, as per Fatunmbi and Salawu [1]. However, due to the diverse characteristics exhibited by non-Newtonian liquids, several empirical rheology formulations have been established based on elastic and nonelastic properties of the fluid, as per Kareem et al. [2]. Oldroyd [3] presented rate-type liquid formulations, Truesdell and Noll [4] considered stress tensor constitutive relations for higher derivatives and velocity gradient symmetric while Rivlin and Ericksen [5] classified viscoelastic liquids based on methodology. Furthermore, in a book, Rajagopal [6] discussed the integral, rate-type and differential type of non-Newtonian fluids. From various non-Newtonian liquids models is the Eyring–Powell fluid, the fluid which describes the viscosity relaxation theory. The liquid is derived fundamentally from fluids kinetic theory and the overall shear stress is based on the combination of type one and two bonds. Therefore, Animasaun et al. [7] considered in a generalized moving plate, a golden thinning molecule and Eyring–Powell fluid motion under buoyancy influence. It was deduced that maximum fluid bulk homogenous species coefficient at the plate occurs at high values of fractional volume and slip velocity parameters. Rahimi et al. [8] employed a collocation technique to give a solution to the Eyring–Powell viscoelastic flow over a linear motioning plate. The study showed that the flow rate is inspired by raising the fluid material term but reduces by enhancing the liquid material term. Salawu et al. [9,10] reported on the variable conductivity of Eyring–Powell hydromagnetic diffusion-reaction fluid flow in a channel. A semi-implicit numerical integration algorithm is adopted for the unsteady state flow, and it can be observed that reactive species solution finite time blowup can be prevented if heat source parameters are consciously managed.

The essential feature of Eyring–Powell fluid alongside with its viscoelastic property has made it valuable to the manufacturing and technological sectors. Coupled with nanofluids, Eyring–Powell fluid usefulness can be enhanced. Nanofluid describes a synthesized liquid that is made up of very small solid suspended nanoparticles of metals with solid fractional volume. Nanofluids exhibits a much better thermal characteristics when compared with traditional fluids, thus, the composition of nanomaterial and base fluids offers a higher thermal conductivity which aids heat transfer, as per Bahiraei et al. [11] and Dadsetani et al. [12]. Bagherzadeh et al. [13] reported that the nanostructural nature of the nanoliquid enhances the chemical, physical, thermal and mechanical properties with an improved convective heat transport and heat conductivity of the fluid. Goshayesh et al. [14] experimentally investigated the two-phase motion and heat transfer of nanofluids in the oscillating heat pipes using the magnetic nanoparticles with kerosene and the iron oxide. Dawar et al. [15] reported the significance of enhancing the thermal characteristics of liquids in various manufacturing and engineering processes especially in cooling and heating processes. In this regard, Khan et al. [16] examined gyrotactic microorganisms in a Walter-B nanomaterial with mixed nonlinear convection and heat diffusion. Thermal diffusion occurs in a vertical unbounded plate, numerical solution taken and the result shows an enhanced temperature field with a rising thermophoresis dispersion and Brownian movement.

Riga plate describes a spanwise aligned array of alternating electrodes and constant magnets fixed on a plane. The initial study of this device was performed by Gailitis and Lielausis [17]. This device produces a crossed electric and magnetic fields capable of stimulating wall parallel Lorentz force such that the fluid flow can be controlled, as per Ahmad et al. [18]. The crucial significance of this device is the reduction in the viscous and pressure drag of submarines as well as prevention of t6he boundary layer separation, Hayat et al. [19]. Hence, the motion of a reactive Eyring–Powell nanoliquid over a an electromagnetic actuator in a non-porous medium was analyzed by Rasool and Zhang [20] while Fatunmbi and Adeosun [21] considered such a problem focusing on the nonlinear radiative flux with heat-mass transfer characteristics and exponential varying viscosity. Despite the applications of nanofluid in improving the characteristics of heat transfer, energy is lost to the ambient in a thermodynamic configuration owing to species reaction, fluid viscosity, mass diffusion and friction force that leads to entropy generation [22]. Due to the existence of the temperature gradient, irreversibility takes place in thermal processes that affect the performance of various thermal devices. Hence, optimization of the engineering devices becomes pertinent which is appropriately carried out using thermodynamic second law.

The thermodynamic second law offers a crucial role in the design and processes of engineering operations, for instance in heat pumps, fire engines, air conditioners, steam power plants and refrigeration works, Pal [22]. As a result, numerous researchers and scientists have analyzed entropy production optimization in Newtonian/non-Newtonian fluids transport subject to various configurations/geometries, boundary conditions and methods. For instance, Muhammad et al. [23] numerically discussed such a concept on a curved surface using the Darcy-Forchheimer model with activation energy. An examination of reactive flow of second-grade nanofluid via the Buongiorno model between two infinite plates in the presence of viscous dissipation and Ohmic heating was scrutinized by Khan et al. [24] using HAM. Goodarzi et al. [25] numerically investigated entropy generation analysis in a mixed convection flow and heat transfer of water considering both laminar and turbulent situations in a square enclosure with the impact of thermal radiation.Similarly, Nakhchi and Esfahani [26] applied the finite volume method to analyze the generation of entropy for the transport of Cu-water nanofluid in a heat exchanger tube having perforated conical rings with large Reynolds number. The analysis showed that Bejan number is a decreasing function of the Reynolds number. Recently, Khosravi et al. [27] engaged the second law of thermodynamics to scrutinize entropy generation for the flow of water-based hybrid nanoliquid passing a wavy cylindrical micro-channel by means of neural networks. Such analysis was also investigated by Pordanjani [28] on free convective flow nanoliquid configured inside a cavity with magnetic field and radiation impact while Yusuf et al. [29] examined such a concept in a hydromagnetic motion of reactive Williamson nanofluid past a convectively heated inclined porous stretching surface with thermal radiation. Fatunmbi and Salawu [30] discussed numerically the entropy production for a magneto-micropolar fluid induced by a nonlinearly stretched sheet and non-uniform magnetic field impact. The analysis showed that the entropy production resulting from heat transfer is high with respect to the radiation parameter.

Various researchers have investigated the flow of Eyring–Powell nanofluids considering diverse parameters of interest on different geometries, assumptions and wall conditions. Nevertheless, the motion of Eyring–Powell nanofluids on magnetized Riga plate has received less attention in the literature. In the available literature, Rasool and Zhang [20] discussed Eyring–Powell nanofluids flow over an impermeable linearly stretching Riga plate with linearized thermal radiation effect. Chu et al. [31] in addition to the parameters investigated by [20] extended such a study by incorporating the effects of nonlinear convection, velocity slip condition and activation energy in the flow field. However, the case of nonlinear stretching Riga plate in a Darcy–Forchheimer porous medium with nonlinear thermal radiation, viscous dissipation coupled with entropy analysis has not been investigated in the the literature. Hence, the current study aims to fill such a gap.

The core focus of the current investigation is to analyze irreversibility for a reactive Eyring–Powell nanoliquid transport over a vertical nonlinearly stretching porous Riga plate with nonlinear thermal radiation. Various reports on nanoliquid Powell–Eyring flow and its significant applications in engineering and manufacturing sectors motivated the study. Complete entropy optimization is considered for four different irreversibilities through thermodynamic second law. The analytical solutions of the problem are carried out using Homotopy analysis method, and compared with Galerkin weighted residual integration technique as presented in tables. Different flow characteristics are examined and the results are graphically presented and discussed appropriately.

2. Problem Modelling

A two-dimensional, incompressible, viscous and steady flow of a non-Newtonian Eyring–Powell nanoliquid past a Riga plate is considered with ($x,y$) coordinate and ($u,v$) being the respective components of velocity (Figure 1). The arrangement of the Riga plate consists of magnets and electrodes which are mounted on the plane as displayed in Figure 1a. The stretching velocity of the plate is indicated by $u=U_w=b{x}^n$ (as depicted in Figure 1b) where n indicates nonlinear stretching term. The combination of the magnetic and electric field of the Riga plate generates the Lorentz force parallel to the wall. The thermal field also features nonlinear radiative flux, viscous dissipation whereas chemical reaction effect is incorporated in the concentration field. The fluid properties are isotropic and constant apart from the thermal conductivity which relies on the temperature. Furthermore, Oberbeck–Boussinesq approximation approach is applied to model the density variation in the body force term of momentum equation.

Following previous authors [32,33], Equation (1) below describes the stress tensor Eyring-Powell liquid.

$$\mathbf{T}=-p\mathbf{l}+\varnothing .$$

(1)

where ${\tau }_{ij}$, the extra stress tensor is written in Equation (2) (see [34])

$${\tau }_{ij}=\frac{1}{\beta }{sinh}^{-1}\left(\frac{1}{\gamma }\frac{\partial u_i}{\partial x_j}\right)+\mu \frac{\partial u_i}{\partial x_j}.$$

(2)

In Equation (2), $\beta$, $\gamma$ and $\mu$ symbolizes material fluid constants and fluid viscosity in that order while the identity tensor is indicated by I. Assuming the second order approximation function is of the form [34].

$${sinh}^{-1}\left(\frac{1}{\gamma }\frac{\partial u_i}{\partial x_j}\right)\cong \frac{1}{\gamma }\frac{\partial u_i}{\partial x_j}-\frac{1}{6}{\left(\frac{1}{\gamma }\frac{\partial u_i}{\partial x_j}\right)}^3,\left|\frac{1}{\gamma }\frac{\partial u_i}{\partial x_j}\right|\ll 1.$$

(3)

In view of the aforementioned assumptions as well as the boundary layer approximation, the governing equations of continuity, momentum, energy and nanoparticle concentration equation are specified below (see [20,32,33,34,35]).

2.1. Governing Equations

Continuity Equation

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

(4)

Momentum Equation

$$\begin{array}{cc}\hfill & \left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=\frac{1}{{\rho }_f}\left({\mu }_f+\frac{1}{\beta \gamma }\right)\frac{{\partial }^{2}u}{\partial y^2}-\frac{1}{2{\rho }_f\beta {\gamma }^3}\frac{{\partial }^{2}u}{\partial y^2}{\left(\frac{\partial u}{\partial y}\right)}^2+\frac{\pi j_0{M}^{★}exp(-\frac{\pi y}{s})}{8{\rho }_f}-\hfill \\ \hfill & \frac{{\mu }_f}{{\rho }_f{K}^{★}}u-F^{★}{u}^2+\frac{1}{{\rho }_f}\left[\left(1-C_{\infty }\right){\rho }_f\left(T-T_{\infty }\right){\beta }_T-\left({\rho }_p-{\rho }_f\right)\left(C-C_{\infty }\right)\right]g,\hfill \end{array}$$

(5)

Energy Equation

$$\begin{array}{cc}\hfill & {\left(\rho c_p\right)}_f\left(u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}\right)=\frac{\partial }{\partial y}\left(k(T)\frac{\partial T}{\partial y}\right)+{\left(\rho c_p\right)}_p\left[\frac{D_T}{T_{\infty }}{\left(\frac{\partial T}{\partial y}\right)}^2+D_B\left(\frac{\partial T}{\partial y}\frac{\partial C}{\partial y}\right)\right]+\hfill \\ \hfill & \left[\left({\mu }_f+\frac{1}{\beta \gamma }\right){\left(\frac{\partial u}{\partial y}\right)}^2-\right.\left.\frac{1}{6\beta {\gamma }^3}{\left(\frac{\partial u}{\partial y}\right)}^4\right]+\frac{{\mu }_f}{c_p{K}^{★}}{u}^2+\frac{{\left(\rho c_p\right)}_f{F}^{★}}{c_p}{u}^3+\frac{16a^{★}}{3b^{★}}\frac{\partial }{\partial y}\left(T^3\frac{\partial T}{\partial y}\right),\hfill \end{array}$$

(6)

Nanoparticles Concentration Equation

$$u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D_B\frac{{\partial }^{2}C}{\partial y^2}+\frac{D_T}{T_{\infty }}\left(\frac{{\partial }^{2}T}{\partial y^2}\right)-k_1\left(C-C_{\infty }\right).$$

(7)

The associated boundary conditions are listed as [18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]

$$\begin{array}{cc}\hfill & u=U_w=b{x}^n,v=V_0,-k_{\infty }\frac{\partial T}{\partial y}=h_T\left(T_f-T\right),-D_B\frac{\partial C}{\partial y}=h_m\left(C_{np}-C\right)aty=0,\hfill \\ \hfill & u\to 0,T\to T_{\infty },C\to C_{\infty },asy\to \infty .\hfill \end{array}$$

(8)

where $h_T=a{x}^{(n-1)/2}$ describes heat transfer coefficient while a is a constant, also, $V_0=V_1{x}^{(n-1)/2}$ represents suction/injection with $V_1$ being a constant [36,37]. More so, $F^{★}=\frac{C_b}{x\sqrt{K_{★}}}$ [38,39] indicates the non-uniform porous medium inertial coefficient while $C_b$ symbolizes the coefficient of drag. Equation (9) is the expression modelling the temperature-dependent thermal conductivity (see [36,37]).

$$k(T)=\frac{k_{\infty }}{T_f-T_{\infty }}\left[(T_f-T_{\infty })+\delta (T-T_{\infty })\right],$$

(9)

here, $k_{\infty }$ denotes upstream thermal conductivity whereas $\delta$ indicates thermal conductivity parameter.

With the following non-dimensional variables, the main equations are dimensionless (see [36,37,38,39,40]).

$$\begin{array}{cc}\hfill & \eta =\sqrt{\frac{b(n+1)x^{n-1}}{2{\nu }_f}}y,\psi =\sqrt{\frac{2{\nu }_{f}b{x}^{n+1}}{n+1}}f(\eta ),\theta (\eta )=\frac{T-T_{\infty }}{T_f-T_{\infty }}\hfill \\ \hfill & \varphi (\eta )=\frac{C-C_{\infty }}{C_{np}-C_{\infty }},u=\frac{\partial \psi }{\partial y},v=-\frac{\partial \psi }{\partial x}.\hfill \end{array}$$

(10)

Substituting Equation (10) into Equations (4)–(7) makes Equation (4) to be satisfied while Equations (5)–(7) in view of Equation (9) simply result to the underlisted equations:

$$\begin{array}{cc}\hfill & \left(1+\alpha \right)f^{‴}-\frac{n+1}{2}\alpha \lambda f^{{}^{″}2}{f}^{‴}+f{f}^{″}-\frac{2n}{n+1}{f}^{{}^{\prime }2}+\frac{2}{n+1}H{e}^{(-B\eta )}+\hfill \\ \hfill & \frac{2}{n+1}{\lambda }_1\left(\theta -R\varphi \right)-\frac{2}{n+1}\left(Da{f}^{\prime }+Fs{f}^{{}^{\prime }2}\right)=0,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & \frac{1}{Pr}\left[1+\delta \theta +Nr{\left(1+({\theta }_b-1)\theta \right)}^3\right]{\theta }^{″}+\frac{3}{Pr}\left[Nr\left({\theta }_b-1\right){\left(1+({\theta }_b-1)\theta \right)}^2\right]{\theta }^{{}^{\prime }2}+\delta {\theta }^{{}^{\prime }2}+f{\theta }^{\prime }+Nt{\theta }^{{}^{\prime }2}+\hfill \\ \hfill & Nb{\theta }^{\prime }{\varphi }^{\prime }+Ec{f}^{{}^{″}2}\left[\left(1+\alpha \right)-\frac{n+1}{6}\alpha \lambda f^{{}^{\prime }2}\right]+\frac{2}{n+1}\left(Da{f}^{{}^{\prime }2}+Fs{f}^{{}^{\prime }3}\right)Ec.\hfill \end{array}$$

(12)

$$\begin{array}{c}\hfill {\varphi }^{″}+Sc\left(f{\varphi }^{\prime }-\frac{2}{n+1}{\gamma }_1\varphi \right)+\frac{Nt}{Nb}{\theta }^{″}.\end{array}$$

(13)

The simplified conditions at the boundary are:

$$\begin{array}{cc}\hfill & f^{\prime }(0)=1,f(0)=F_w,\theta (0)=1+\frac{1}{B_1}{\theta }^{\prime }(0),\varphi (0)=1+\frac{1}{B_2}{\varphi }^{\prime }(0),\hfill \\ \hfill & f^{\prime }(\infty )=0,\theta (\infty )=0,\varphi (\infty )=0.\hfill \end{array}$$

(14)

2.2. The Entropy Generation Equation

By the application of thermodynamics second law, the volumetric rate of entropy production for the Eyring–Powell fluid influenced by viscous dissipation, nonlinear thermal radiation and diffusion effect is modelled as [41,42,43].

$$S_{Gn}=S_{ht}+S_v+S_d,$$

(15)

where $S_{Gn}$ is the volumetric entropy generation, $S_{ht}$ typifies entropy production owing to heat transfer or the conduction effect, $S_v$ represents viscous dissipation entropy production whereas the last term indicates the entropy production with mass transfer effects across the finite difference in temperature and concentration. Thus, $S_{Gn}$ can be written as:

$$\begin{array}{cc}\hfill & S_{Gn}=\frac{1}{T^2}\left(k+\frac{16a^{★}{T}^3}{3b^{★}k}\right){\left(\frac{\partial T}{\partial y}\right)}^2+\frac{1}{{\rho }_f{c}_p}\left[\left({\mu }_f+\frac{1}{\beta \gamma }\right){\left(\frac{\partial u}{\partial y}\right)}^2-\right.\hfill \\ \hfill & \left.\frac{1}{6\beta {\gamma }^3}{\left(\frac{\partial u}{\partial y}\right)}^4\right]+\frac{{\mu }_f}{T{K}^{★}}{u}^2+\left[\frac{R{D}_B}{C}{\left(\frac{\partial C}{\partial y}\right)}^2+\frac{R{D}_B}{T}\left(\frac{\partial C}{\partial y}\frac{\partial T}{\partial y}\right)\right].\hfill \end{array}$$

(16)

where $S_G^{‴}$ indicates characteristic entropy generation which is expressed as

$$S_G^{‴}=k_{\infty }\frac{{(T_w-T_{\infty })}^2(n+1)}{2T_{\infty }^2{x}^2}$$

(17)

Similarly, the dimensionless entropy generation is expressed as:

$$\begin{array}{cc}\hfill & N_{Gs}=\frac{Re\left[1+\delta \theta +Nr{(1+\left({\theta }_b-1\right)\theta )}^3\right]}{{(1+\left({\theta }_b-1\right)\theta )}^2}{\theta }^{{}^{\prime }2}+\frac{ReBr}{\left(1+({\theta }_b-1)\theta \right)({\theta }_b-1)}\left[(1+\Omega )f^{{}^{″}2}-\left(\frac{\lambda \alpha (n+1)}{6}\right)f^{{}^{″}4}\right]+\hfill \\ \hfill & \frac{ReBrDa}{\left(1+({\theta }_b-1)\theta \right)({\theta }_b-1)}{f}^{{}^{\prime }2}+\frac{Re{\gamma }_2({\varphi }_b-1)}{({\theta }_b-1)}\left[\frac{({\varphi }_b-1){\varphi }^{{}^{\prime }2}}{\left(1+({\varphi }_b-1)\varphi \right)({\theta }_b-1)}+\frac{{\theta }^{\prime }{\varphi }^{\prime }}{1+({\theta }_b-1)\theta }\right]=0.\hfill \end{array}$$

(18)

where the term $N_{Gs}=\frac{S{G}_n}{S_G^{‴}}\left(=N_h+N_v+N_m\right)$ showcases the entropy production number. In this expression, $N_h$ indicates the first term in the RHS of Equation (18) which implies entropy production owing to heat transfer while $N_v$ is sum of the second and third terms in the RHS of Equation (18) and describes entropy generation due to frictional heating, whereas the last term $N_m$ implies entropy generation as a result of mass transfer. Similarly, the Bejan number $Be$ is also of paramount interest to the engineers, it defines the ratio of the heat and mass transfer irreversibility to that of entropy generation number. Mathematically, $Be$ can be expressed as

$$Be=\frac{N_h+N_m}{N_{Gs}}=\frac{1}{1+\mathrm{Y}},$$

(19)

where $\mathrm{Y}=\frac{N_v}{(N_h+N_m)}$ which shows the irreversibility ratio. The dominance of frictional irreversibility over the heat and mass transfer occurs when $\mathrm{Y}>1$, the converse is true when $0\le \mathrm{Y}<1$ whereas equal contribution of frictional heating irreversibility to that of heat and mass transfer ($\mathrm{Y}=1$). The Bejan number ranges in the interval $0\le Be\le 1$. When $Be=1$, heat transfer and mass irreversibility overrides the frictional heating effect whereas the opposite occurs when $Be=0$. The situation where $Be=0.5$ indicates equal contributions.

The emerging parameters from the governing equations are described as:

$$\begin{array}{cc}\hfill & {\theta }_b=\frac{T_f}{T_{\infty }},\alpha =\frac{1}{{\mu }_f\beta \gamma },Da=\frac{{\nu }_f}{b{K}^{★}},H=\frac{\pi j_0{M}^{★}}{8{\rho }_f{b}^2{x}^{2n-1}(n+1)},Pr=\frac{{\mu }_f{c}_p}{k_{\infty }},Fs=\frac{Cb}{\sqrt{K^{★}}},Nr=\frac{16a^{★}{T}_{\infty }^3}{3b^{★}{k}_{\infty }},\hfill \\ \hfill & G{r}_x=\frac{g{\beta }_T\left(1-C_{\infty }\right)\left(T_f-T_{\infty }\right)x^3}{{\nu }_f^2},\lambda =\frac{b^3{x}^{3n-1}}{4{\nu }_f{\gamma }^2},Ec=\frac{U_w^2}{C_p(T_w-T_{\infty })},Sc=\frac{{\nu }_f}{Dm},\hfill \\ \hfill & Nb=\frac{{(\rho c_p)}_p{D}_B\left(C_w-C_{\infty }\right)}{{\nu }_f{(\rho c_p)}_f},R{e}_x=\frac{c{x}^{n+1}}{{\nu }_f},{\lambda }_1=\frac{g{\beta }_T\left(1-C_{\infty }\right){\rho }_{f\infty }\left(T_f-T_{\infty }\right)}{{\rho }_f{c}^2{x}^{2n-1}}=\frac{G{r}_x}{R{e}_x^2},\hfill \\ \hfill & B_1=\frac{a}{k_{\infty }}\sqrt{\frac{2{\nu }_f}{b(n+1)}},B_2=\frac{a}{D_B}\sqrt{\frac{2{\nu }_f}{b(n+1)}},Fw=\sqrt{\frac{2V_1^2}{b{\nu }_f(n+1)}},B=\frac{2\pi }{s}\sqrt{\frac{2{\nu }_f}{b(n+1)x^{n+1}}},\hfill \\ \hfill & R=\frac{\left({\rho }_p-{\rho }_f\right)\left(C_w-C_{\infty }\right)}{{\beta }_T\left(1-C_{\infty }\right)\left(T_f-T_{\infty }\right){\rho }_f},{\gamma }_1=\frac{k_1}{b{x}^{n-1}},Nt=\frac{{(\rho c_p)}_p{D}_T\left(T_f-T_{\infty }\right)}{{\nu }_f{T}_{\infty }{(\rho c_p)}_f}.\hfill \end{array}$$

(20)

The particular quantities of interest to the engineers are the skin friction coefficient $c_{fx}$, the local Nusselt number $N_{ux}$ and the local Sherwood number $S_{hx}$ which are orderly described as

$$c_{fx}=\frac{{\tau }_w}{{\rho }_f{U}_w^2},N_{ux}=\frac{x{q}_w}{k_{\infty }\left(T_f-T_{\infty }\right)},S_{hx}=\frac{x{q}_m}{D_B\left(C_w-C_{\infty }\right)},$$

(21)

in Equation (21), ${\tau }_w$ defines the shear stress, $q_w$($q_m$) defines heat (mass) flux at the surface. These are sequentially expressed as

$${\tau }_w={\left.\left[\left({\nu }_f+\frac{1}{\rho \beta \gamma }\right)\frac{\partial u}{\partial y}+\frac{1}{6\beta {\gamma }^3}{\left(\frac{\partial u}{\partial y}\right)}^3\right]\right|}_{y=0},q_w=-{\left.\left(k+\frac{16T^3\sigma }{3k^{★}}\right)\frac{\partial T}{\partial y}\right|}_{y=0},q_m=-D_B{\left.\frac{\partial C}{\partial y}\right|}_{y=0},$$

(22)

invoking Equations (10) and (22) in (21) results to Equations (23)–(25)

$$c_{fx}={\left(\frac{n+1}{2}\right)}^{1/2}\left[(1+\alpha )f^{″}(0)-\frac{1}{3}\alpha \lambda f^{″}{(0)}^3\right]R{e}_x^{-1/2},$$

(23)

$$N_{ux}=-{\left(\frac{n+1}{2}\right)}^{1/2}\left[1+Nr{\left(1+({\theta }_b-1)\theta (0)\right)}^3\right]R{e}_x^{1/2}{\theta }^{\prime }(0),$$

(24)

$$S_{hx}=-{\left(\frac{n+1}{2}\right)}^{1/2}R{e}_x^{\frac{1}{2}}{\varphi }^{\prime }(0).$$

(25)

3. Method of Solution

The boundary value problem (11)–(14) have been analytically tacked via the Homotopy Analysis Method (HAM). HAM is an analytic approximation recursive technique which yields a series sum solution and offers a simple way to guarantee the convergence of the solution series. It also provides freedom to select auxiliary linear operator, initial guess and solution expression for higher order. For details concerning the method see [44,45]. For the solutions of Equations (11)–(13), the following initial guesses and linear operators are considered

$$f_0=1+F_w-e^{-\eta },{\theta }_0=\frac{B_1{e}^{-\eta }}{1+B_1},{\varphi }_0=\frac{B_2{e}^{-\eta }}{1+B_2},$$

(26)

$$L_f=\frac{d^3}{d{\eta }^3}+\frac{d^2}{d{\eta }^2},L_{\theta }=\frac{d^2}{d\eta }+\frac{d}{d\eta },L_{\varphi }=\frac{d^2}{d\eta }+\frac{d}{d\eta },$$

(27)

with the properties that

$$L_f(D_1+D_2\eta +D_3{e}^{-\eta })=0,L_{\theta }(D_4+D5e^{-\eta })=0,and{L}_{\varphi }(D_6+D_7{e}^{-\eta })=0.$$

(28)

where $D_j(j=1,2,\cdots ,7)$ are unknown constants to be obtained from the wall conditions. For the embedding parameter $q\in [0,1]$, one can express the zeroth order deformation problem as

$$(1-q)L_f\left[\widehat{f}(\eta ;q)-f_0(\eta )\right]=q{\hslash }_f{\aleph }_f\left[\widehat{f}(\eta ;q),\widehat{\theta }(\eta ;q),\widehat{\varphi }(\eta ;q)\right]$$

(29)

$$(1-q)L_{\theta }\left[\widehat{\theta }(\eta ;q)-{\theta }_0(\eta )\right]=q{\hslash }_{\theta }{\aleph }_{\theta }\left[\widehat{f}(\eta ;p),\widehat{\theta }(\eta ;q),\widehat{\varphi }(\eta ;q)\right]$$

(30)

$$(1-q)L_{\varphi }\left[\widehat{\varphi }(\eta ;q)-{\varphi }_0(\eta )\right]=q{\hslash }_{\varphi }{\aleph }_{\varphi }\left[\widehat{f}(\eta ;q),\widehat{\theta }(\eta ;q),\widehat{\varphi }(\eta ;q)\right]$$

(31)

The non-zero auxiliary parameters are described as ${\hslash }_f,{\hslash }_{\theta }$ and ${\hslash }_{\varphi }$ and ${\aleph }_f,{\aleph }_{\theta }$, and ${\aleph }_{\varphi }$ are nonlinear operators. When $q=0$ and $q=1$ we have

$$\widehat{f}(\eta ;0)=f_0(\eta ),\widehat{f}(\eta ;1)=f(\eta ),\widehat{\theta }(\eta ;0)={\theta }_0(\eta ),\widehat{\theta }(\eta ;1)=\theta (\eta ),\widehat{\varphi }(\eta ;0)={\varphi }_0(y),\widehat{\varphi }(\eta ;1)=\varphi (\eta )$$

(32)

The variation of q from 0 to 1 implies that $\widehat{f}(\eta ;q),\widehat{\theta }(\eta ;q)$ and $\widehat{\varphi }(\eta ;q)$ vary from $f_0(\eta ),{\theta }_0(\eta )$ and ${\varphi }_0(\eta )$ to $f(\eta ),\theta (\eta )$ and $\varphi (\eta )$. Now, the Taylor series expansion of $\widehat{f}(\eta ;q),\widehat{\theta }(\eta ;q)$ and $\widehat{\varphi }(\eta ;q)$ with respect to q yield

$$f(\eta ;q)=f_0(\eta )+\sum _{m=1}^{\infty }{f}_m(\eta )q^m$$

(33)

$$\theta (\eta ;q)={\theta }_0(\eta )+\sum _{m=1}^{\infty }{\theta }_m(\eta )q^m$$

(34)

$$\varphi (\eta ;q)={\varphi }_0(\eta )+\sum _{m=1}^{\infty }{\varphi }_m(\eta )q^m$$

(35)

where

$$f_m(\eta )=\frac{1}{m!}{\left.\frac{{\partial }^{m}f(\eta ;q)}{\partial q^m}\right|}_{q=0},{\theta }_m(\eta )=\frac{1}{m!}{\left.\frac{{\partial }^m\theta (\eta ;q)}{\partial q^m}\right|}_{q=0},{\varphi }_m(\eta )=\frac{1}{m!}{\left.\frac{{\partial }^m\varphi (\eta ;q)}{\partial q^m}\right|}_{q=0}$$

(36)

Assuming that auxiliary parameters, initial guesses, auxiliary linear operators are properly chosen such that Equations (33)–(35) converge at $q=1$, we obtain

$$f(\eta )=f_0(\eta )+\sum _{m=1}^{\infty }{f}_m(\eta )$$

(37)

$$\theta (\eta )={\theta }_0(\eta )+\sum _{m=1}^{\infty }{\theta }_m(\eta )$$

(38)

$$\varphi (\eta )={\varphi }_0(\eta )+\sum _{m=1}^{\infty }{\varphi }_m(\eta )$$

(39)

The corresponding mth order deformation problems are:

$$L_f\left[f_m(\eta )-X_m{f}_{m-1}(\eta )\right]={\hslash }_f{R}_{f,m-1}(\eta )$$

(40)

$$L_{\theta }\left[{\theta }_m(\eta )-X_m{\theta }_{m-1}(\eta )\right]={\hslash }_{\theta }{R}_{\theta ,m-1}(\eta )$$

(41)

$$L_{\varphi }\left[h_m(\eta )-X_m{h}_{m-1}(\eta )\right]={\hslash }_{\varphi }{R}_{\varphi ,m-1}(\eta )$$

(42)

subjected to homogeneous boundary conditions

$$f_m(0)=f_m^{\prime }(0)=f_m^{\prime }(\infty )={\theta }_m^{\prime }(0)-B_1{\theta }_m(0)={\theta }_m(\infty )={\varphi }_m(0)={\varphi }_m(\infty )=0$$

(43)

where

$$R_{f,m-1}(\eta )=\frac{1}{(m-1)!}{\left.\frac{{\partial }^{(m-1)}{\aleph }_f[f(\eta ;q),\theta (\eta ;q),\varphi (\eta ;q)]}{\partial q^{(m-1)}}\right|}_{q=0}$$

(44)

,

$$R_{\theta ,m-1}(\eta )=\frac{1}{(m-1)!}{\left.\frac{{\partial }^{(m-1)}{\aleph }_{\theta }[f(\eta ;q),\theta (\eta ;q),\varphi (\eta ;q)]}{\partial q^{(m-1)}}\right|}_{q=0},$$

(45)

$$R_{\varphi ,m-1}(\eta )=\frac{1}{(m-1)!}{\left.\frac{{\partial }^{(m-1)}{\aleph }_{\varphi }[f(\eta ;q),\theta (\eta ;q),\varphi (\eta ;q)]}{\partial q^{(m-1)}}\right|}_{q=0}.$$

(46)

$X_m=0$ for $m\le 1$ and $X_m=1$ for $m>1$. The general solutions can be written as:

$$f_m(\eta )=f_m^{*}(\eta )+D_1+D_2\eta +D_3{e}^{-\eta },{\theta }_m(\eta )={\theta }_m^{*}(\eta )+D_4+D5e^{-\eta },{\varphi }_m(\eta )={\varphi }_m^{*}(\eta )+D_6+D_7{e}^{-\eta }$$

(47)

$f_m^{*}(\eta ),{\theta }_m^{*}(\eta ),{\varphi }_m^{*}(\eta )$ are the particular solutions and $D_j(j=1,2,\cdots ,7)$ are obtained from the conditions in Equation (43). Mathematical symbolic package, MATHEMATICA 11.3 is employed to solve Equations (40)–(42) in turns in the order of $m=1,2,\cdots$ (see the code in the algorithm in Supplemantary Materials).

Convergence of HAM

Proper selection of the values of auxiliary parameters contributes to the convergence of HAM. To ascertain the convergence of the obtained results, the interval of ℏ curves for the function $f^{″}(0)$, ${\theta }^{\prime }(0)$ and ${\varphi }^{\prime }(0)$ with default values of $\lambda =n=H={\lambda }_1=R=Da=Fs=\delta =Nt=Ec={\gamma }_1=\alpha =Fw=0.1$, $B_1=0.2$, $Pr=Sc=1,Nb=0.5,Nr=0.3,{\theta }_b=1.5,B=2.0$ at 16th order of approximation is displayed in Figure 2. The admissible range of ${\hslash }_f,{\hslash }_{\theta }$ and ${\hslash }_{\varphi }$ is found at $-0.92\le {\hslash }_f\le -0.14$, $-1.04\le {\hslash }_{\theta }\le -0.05$ and $-1.14\le {\hslash }_{\varphi }\le 0.05$, respectively. For the realization of the optimal values for ${\hslash }_f,{\hslash }_{\theta }$ and ${\hslash }_{\varphi }$, square residual error minimization is adopted as follows:

$${\Delta }_m^f=\frac{1}{K+1}\sum _{j=0}^K{\left[{\aleph }_f{\left(\sum _{k=0}^m{f}_k\left(\eta \right),\sum _{k=0}^m{\theta }_k\left(\eta \right),\sum _{k=0}^m{h}_k\left(\eta \right)\right)}_{\eta =j\Delta x}\right]}^2$$

(48)

$${\Delta }_m^{\theta }=\frac{1}{K+1}\sum _{j=0}^K{\left[{\aleph }_{\theta }{\left(\sum _{k=0}^m{f}_k\left(\eta \right),\sum _{k=0}^m{\theta }_k\left(\eta \right),\sum _{k=0}^m{h}_k\left(\eta \right)\right)}_{\eta =j\Delta x}\right]}^2$$

(49)

$${\Delta }_m^h=\frac{1}{K+1}\sum _{j=0}^K{\left[{\aleph }_{\varphi }{\left(\sum _{k=0}^m{f}_k\left(\eta \right),\sum _{k=0}^m{\theta }_k\left(\eta \right),\sum _{k=0}^m{h}_k\left(\eta \right)\right)}_{\eta =j\Delta x}\right]}^2$$

(50)

Following Liao [46],

$${\Delta }_m^t={\Delta }_m^f+{\Delta }_m^{\theta }+{\Delta }_m^{\varphi }$$

(51)

where $\Delta x=\frac{10}{K},K=50$ and ${\Delta }_m^t$ is the total square residual error. At the fourth order of approximation, ${\Delta }_m^t$ is minimized to obtain optimal values of ${\hslash }_f,{\hslash }_{\theta }$ and ${\hslash }_{\varphi }$ as follows

$$\frac{\partial {\Delta }_m^t}{\partial {\hslash }_f}=0,\frac{\partial {\Delta }_m^t}{\partial {\hslash }_{\theta }}=0,\frac{\partial {\Delta }_m^t}{\partial {\hslash }_{\varphi }}=0$$

(52)

Equation (52) is solved simultaneously and the values of auxiliary parameters are obtained as ${\hslash }_f=-0.74597$, ${\hslash }_{\theta }=-0.77980$ and ${\hslash }_h=-0.57914$. The convergence of solution approximated at diverse orders is recorded in

Table 1. HAM solution Convergence approximated at diverse orders when ${\hslash }_f=-0.74597$, ${\hslash }_{\theta }=-0.77980$ and ${\hslash }_{\varphi }=-0.57914$.

m	$-{\mathit{f}}^{″}(0)$	$-{\mathit{\theta }}^{\prime }(0)$	$-{\mathit{\varphi }}^{\prime }(0)$
2	0.75237	0.12555	0.15604
4	0.73895	0.12408	0.15635
6	0.74000	0.12553	0.15648
8	0.74057	0.12548	0.15635
10	0.74055	0.12535	0.15633
12	0.74051	0.12537	0.15634
14	0.74051	0.12538	0.15634
16	0.74051	0.12538	0.15634
18	0.74051	0.12538	0.15634
20	0.74051	0.12538	0.15634

. In

Table 2. Validation of Homotopy Analysis Method (HAM) results at sixteenth order with results obtained via Galerkin weighted residual method (GWRM).

Parameters	Values	$-{\mathit{\theta }}^{\prime }(0)$		$-{\mathit{\varphi }}^{\prime }(0)$
		HAM	GWRM	HAM	GWRM
$B_1$	0.1	0.07491	0.07491	0.15716	0.15716
	0.3	0.16054	0.16053	0.15582	0.15582
	0.5	0.20488	0.20485	0.15525	0.15524
${\lambda }_1$	0.1	0.12538	0.12538	0.15634	0.15634
	0.3	0.12796	0.12794	0.15682	0.15682
	0.5	0.12983	0.12981	0.15721	0.15721
R	0.1	0.12538	0.12538	0.15634	0.15634
	0.3	0.12525	0.12525	0.15632	0.15632
	0.5	0.12512	0.12511	0.15629	0.15629
$\alpha$	0.1	0.12538	0.12538	0.15634	0.15634
	0.3	0.12601	0.12601	0.15679	0.15679
	0.5	0.12646	0.12646	0.15717	0.15717
$\lambda$	0.1	0.12538	0.12538	0.15634	0.15634
	0.3	0.12538	0.12537	0.15634	0.15634
	0.5	0.12537	0.12537	0.15633	0.15633
H	0.1	0.12538	0.12538	0.15634	0.15634
	0.3	0.12692	0.12692	0.15666	0.15666
	0.5	0.12828	0.12828	0.15698	0.15698

, the obtained results of ${\theta }^{\prime }(0),{\varphi }^{\prime }(0)$ using HAM is validated via Galerkin weighted residual method (GWRM) and a good agreement is found.

4. Results and Discussion

For better understanding of the impact of the embedded parameters on the dimensionless quantities, several graphs have been constructed in this section to showcase such reactions. Figure 3 depicts the reaction of the fluid parameter $\alpha$ on the velocity field in the presence of the buoyancy ratio term R. It is obvious that an uplift in both $\alpha$ and R heighten the velocity field and accelerates the fluid motion. The quantity $\alpha$ varies inversely to the fluid viscosity, hence, a rise $\alpha$ implies a reduction of viscosity and at such, the fluid motion enhances.

In a similar situation, the plot of describing the velocity field against $\eta$ for diverse values of H (modified Hartmann number) in the existence of the dimensionless parameter B reveals that the velocity field accelerates with a hike in H as shown in Figure 4. Actually, this pattern agrees well with the physical principles of the problem in the sense that $H>0$ signifies assisting flow mechanism on the field of velocity. The generated Lorentz forces by the Riga plate parallel to the surface promote surface tension which enhances the fluid motion. However, uplifting the dimensionless parameter B acts contrarily due to a rise in the fluid viscosity and, as such, a resistance to fluid motion is created. Thus, improvement in B decelerates the velocity profile together with the associated hydrodynamic boundary structure as displayed in Figure 4. The plot in Figure 5 depicts the response of the velocity field with variations in ${\lambda }_1$ (mixed convection term) and $F_w$. Growing values of ${\lambda }_1$ boosts fluid motion due to an enhancement in the buoyancy force whereas a hike in the suction term reduces the motion of the fluid. The fluid motion also drops with an increase in the Darcy term $Da$ as described in Figure 6. The reaction of the surface convection term $B_1$ (Biot number) on the thermal field is sketched in Figure 7. Vividly, there is an expansion in the size of the thermal boundary layer with a lift in $B_1$ and in consequence, the temperature distribution improves as shown in the figure. A boost in $B_1$ strengthens the transport of heat coefficient such that the temperature field escalates. Likewise, there is an enhancement in the temperature field with escalating values $\delta$ and $Nb$ as demonstrated in Figure 8. In a physical sense, a boost in $Nb$ compels a rise in the transport of the nanomaterials and at such, the thermal field expands.

The attribute of the thermal field to variations in $Nr$ and ${\theta }_b$ is constructed in Figure 9. The thermal field increases with enhancement in the radiation parameter. Likewise, a growth in the temperature parameter ${\theta }_b$ favours the operating temperature difference $T_f-T_{\infty }$ such that the surface temperature improves as illustrated in Figure 9. The impact of chemical reaction ${\gamma }_1$ together with Brownian motion $Nb$ effect in respect to concentration field are highlighted in Figure 10. Vividly, the plot reveals that a hike in the value of ${\gamma }_1$ compels the concentration boundary layer to decline significantly. Similarly, a growth in the magnitude of the Brownian motion $Nb$ leads to a reduction in the concentration profile due to the fact that higher Brownian motion creates a stronger collision of the fluid particles.

The plot of the concentration profile versus $\eta$ for variation in $Sc$ in the existence of $Nt$ is described in Figure 11. The plot reveals that the concentration field decreases as $Sc$ rises whereas the converse is the case when $Nt$ escalates. Figure 12 explains the reaction of $N_{Gs}$ (entropy generation number) for changes in the values of ${\theta }_b$ and $\alpha$. Obviously, increasing values of ${\theta }_b$ lowers $N_{Gs}$. This trend agree well with the report discussed by Sithole et al. [47]. On the other hand, the entropy production escalates with higher values of the fluid parameter $\alpha$. Furthermore, a rise in the value of $Ec$ raises the entropy generation due to a hike in the frictional forces as $Ec$ increases (see Figure 13). Likewise, an uplift in the magnitude of the Darcy number ($Da$) strengthens the production of entropy ($N_{Gs}$) as demonstrated in Figure 13.

The entropy generation $N_{Gs}$ also abounds with a raise in $Nr$ (radiation term) as described in Figure 14. An increase in $Nr$ enhances the temperature which in turn affects $N_{Gs}$. Furthermore, a rise in $B_1$ boosts the entropy production $N_{Gs}$. This fact is described in Figure 14. Moreover, Figure 15 demonstrates the influence of $Nr$ and ${\theta }_b$ on the Bejan number $Be$. In this figure, it is noticed that $Be$ improves with an uplift in both $Nr$ and ${\theta }_b$. The implication here is that the entropy production resulting from heat and mass transfer dominate that of the frictional heating effect with rising values of $Nr$ and ${\theta }_b$. On the contrary, Figure 16 elucidates that the frictional heating entropy production overrides that of heat and mass transfer when $Ec$ and $Da$ are raised due to a fall in $Be$.

5. Conclusions

A mathematical model has been developed to analyze entropy generation in a reactive Eyring–Powell nanofluid flow over a vertical Riga plate in a porous device. Incorporated in the model are the influences of nonlinear thermal radiation, thermophoresis, viscous dissipation, Brownian motion coupled with varying thermal conductivity in the presence of convective thermal and concentration wall conditions. An analytical solution via Homotopy analysis method has employed to obtain solutions to the governing transport equations. Furthermore, the obtained solutions via HAM have been validated by the Galerkin weighted residual method and a perfect agreement exists for the selected parameters. The parametric evaluations of the effects of the emerging parameters on the dimensionless quantities are graphically displayed and consequently deliberated in the study. The analysis has revealed that:

• Augmenting the modified Hartmann number, mixed convection and buoyancy terms enlarge the hydrodynamic boundary layer leading to enhancement of the velocity field whereas the hydrodynamic boundary structure shrinks with incremental values of suction and Darcy parameters.
• The thermal field expands with growth in the surface convection term (Biot number), thermal conductivity, radiation, thermophoresis as well as Brownian motion and temperature ratio parameters whereas the concentration profile behaves conversely when Schmidt number, Brownian motion and chemical reaction terms increase.
• Entropy production is high when Eckert number, radiation term as well as Darcy number increases while such a pattern changes with an uplift in the temperature ratio term.
• The consequences of increasing radiation and temperature ratio parameters are an improvement of Bejan number which in turn leads to the dominance of entropy production due to heat and mass transfer over that of frictional heating effect. This trend is, however, reversed with advancing Eckert and Darcy numbers as the duo deplete Bejan number.