Entropy Generation Analysis for Combustible Third Grade Fluid Flow through a Slant Channel

1. Introduction

A good number of researchers have paid renewed interest and attention to the flow of reactive fluids because due to technological advancements. In particular, energy transfer and variable viscosity property of a fluid have been of vast interest among researchers [1,2,3,4,5,6,7,8]. Industrial and technological applications are found in areas such as atmospheric flows, thermal regulation, cooling of electronic devices, nuclear reactors, thermal hydraulics and many others other devices have proliferated the engineering community.

Meanwhile, entropy generation minimization in thermo-fluid flows is critically important in improving efficiency of machines and processes as explained by Bejan [9]. Therefore, efficient energy conversion and management is key in any thermal processes. In the industry, issues like productivity, sustainability and competitiveness require engineering solutions, and such solutions are heavily reliant on mathematical models. Modern economies are being urged to curtail their energy dependence on fossil fuels in order to mitigate continuous depletion of the ozone layer that is causing acceleration of global warming. The deleterious consequences of global warming like runaway fires, rising sea levels and melting of glaciers are occurring at unsustainable rates in the present day. As more clean renewable energy forms are being sought to replace the fossil fuels, it has become more and more urgent for economic systems to conserve more and more energy for continued sustainability. Thus, the need for continued research on entropy generation minimization cannot be over-emphasized.

The assumption of the no-slip boundary condition has been the norm in most studies involving impervious and non-polished surfaces. Many researchers are of the opinion that the no-slip boundary condition is an hypothesis. In fact evidences of slip of a fluid on a solid surface were reported by several authors including Mathews and Hill [10] and Zhu and Granick [11]. Rundora and Makinde [12,13] reported this effect on channel flows saturated with porous barriers. Das et al. [14] reported a magnetized flow with flud slippage. The aim of this present study is to formulate a theoretical model where we investigate the combined effects of Navier slip, variable viscosity, porous medium permeability and convective boundary conditions on entropy generation in a steady flow of a reactive third-grade fluid through an inclined channel filled with a porous saturated medium. A survey of literature revealed that the scope of our study herein has not been fully accounted for in the previous studies.

2. Mathematical Formulation

Consider a steady flow of an incompressible variable viscosity, reactive third grade fluid through an inclined channel filled with a homogeneous and isotropic porous medium as shown in Figure 1 The lower wall of the channel is subjected to convective heat exchange with the surrounding medium, while the upper wall is subjected to a constant heat flux. It is assumed that the convective heat exchange with the ambient follows Newton’s law of cooling. The flow is assumed to be induced by an applied axial constant pressure gradient and bouyancy force. Neglecting the reacting viscous fluid consumption the model equations emanating from the momentum and heat balance can be written as [12,15,16]:

$$0=-{\displaystyle \frac{d\overline{P}}{d\overline{x}}}+{\displaystyle \frac{d}{d\overline{y}}}\left(\overline{\mu }\left(T\right){\displaystyle \frac{d\overline{u}}{d\overline{y}}}\right)+6{\beta }_3{\displaystyle \frac{d^2\overline{u}}{d{\overline{y}}^2}}{\left({\displaystyle \frac{d\overline{u}}{d\overline{y}}}\right)}^2-{\displaystyle \frac{\overline{\mu }\left(T\right)\overline{u}}{K}}+\rho g\beta (T-T_0)sin\left(\omega \right),$$

(1)

$$0=k{\displaystyle \frac{d^{2}T}{d{\overline{y}}^2}}+{\left({\displaystyle \frac{d\overline{u}}{d\overline{y}}}\right)}^2\left[\overline{\mu }\left(T\right)+2{\beta }_3{\left({\displaystyle \frac{d\overline{u}}{d\overline{y}}}\right)}^2\right]+{\displaystyle \frac{\overline{\mu }\left(T\right){\overline{u}}^2}{K}}+Q{C}_{0}A{\left({\displaystyle \frac{hT}{vl}}\right)}^m{e}^{-{\displaystyle \frac{E}{RT}}},$$

(2)

with appropriate boundary conditions

$$\begin{array}{c}\hfill \overline{y}=-a:{\lambda }_1\overline{u}=\overline{\mu }\left(T\right){\displaystyle \frac{d\overline{u}}{d\overline{y}}},-k{\displaystyle \frac{dT}{d\overline{y}}}=h_1[T-T_0],\\ \hfill \overline{y}=a:{\lambda }_2\overline{u}=\overline{\mu }\left(T\right){\displaystyle \frac{d\overline{u}}{d\overline{y}}},-k{\displaystyle \frac{dT}{d\overline{y}}}=q_w.\end{array}$$

(3)

$\overline{P}$ is the modified fluid pressure, $\overline{x}$ and $\overline{y}$ are the axial and normal coordinates to the inclined channel, $\overline{u}$ is the fluid velocity, $h_1$ is the heat transfer coefficient at the lower plate, $T_0$ is the fluid initial temperature, T is the fluid temperature, g is the acceleration due to gravity, $\omega$ is the angle of inclination, ${\beta }_3$ is the third-grade material coefficient, k is the thermal conductivity, K is the Porous medium permeability, $\overline{\mu }$ is the fluid dynamic viscosity, $\rho$ is the density, Q is the heat generated internally, $C_0$ is the initial concentration of the reactant species, A is the reaction rate constant, ${\lambda }_1$ and ${\lambda }_2$ are the slip coefficients at the lower and upper channel walls respectively, $q_w$ is the constant heat flux, E is the activation energy, h is the Boltzmann’s constant, l is the Planck’s number, R is the universal gas constant, v is the vibration frequency, $\beta$ is the volumetric thermal expansion coefficient, m is the numerical exponent such that the three values represent numerical exponents for sensitised, Arrhenius and biomelecular kinetics respectively as $m\in \{-2,0,0.5\}$ [17,18]. The temperature dependent viscosity $\overline{\mu }\left(T\right)$ can be expressed as

$$\overline{\mu }\left(T\right)={\mu }_0{e}^{-b(T-T_0)}$$

(4)

where ${\mu }_0$ is the initial fluid viscosity at temperature $T_0$ and b is the viscosity variation parameter. Introducing the following non-dimensional variables into Equations (1)–(4),

$$\begin{array}{c}\hfill y=\frac{\overline{y}}{a},x=\frac{\overline{x}}{a},f=\frac{\overline{u}\rho a}{{\mu }_0},\alpha =\frac{bR{T}_0^2}{E},G=-\frac{dP}{dx},\Omega =\frac{{\beta }_3{\mu }_0}{{\rho }^2{a}^4},Da=\frac{K}{a^2},\\ \hfill P=\frac{a^2\rho \overline{P}}{{\mu }_0^2},ϵ=\frac{R{T}_0}{E},\theta =\frac{E(T-T_0)}{R{T}_0^2},Gr=\frac{{\rho }^{2}g\beta a^{3}R{T}_0^2}{{\mu }_0^{2}E},Br=\frac{{\mu }_0^{3}E}{a^2{\rho }^{2}kR{T}_0^2},\\ \hfill {\delta }^2=\frac{1}{Da},{\beta }_1=\frac{{\mu }_0}{{\lambda }_{1}a},{\beta }_2=\frac{{\mu }_0}{{\lambda }_{2}a},Bi=\frac{a{h}_1}{k},\lambda =\frac{QEA{a}^2{C}_0}{kR{T}_0^2}{\left(\frac{h{T}_0}{\nu l}\right)}^m{e}^{-\frac{E}{R{T}_0}},\end{array}$$

(5)

we obtain the non-dimensional governing equations

$$e^{-\alpha \theta }{f}^{\prime \prime }+6\Omega f^{\prime \prime }{f}^{\prime 2}-\alpha e^{-\alpha \theta }{f}^{\prime }{\theta }^{\prime }-\delta e^{-\alpha \theta }f+Gr\theta Sin\left(\omega \right)+G=0,$$

(6)

$${\theta }^{\prime \prime }+Br{f}^{\prime 2}(e^{-\alpha \theta }+2\Omega f^{\prime 2})+\delta Br{e}^{-\alpha \theta }{f}^2+\lambda {(1+ϵ\theta )}^{m}exp\left(\frac{\theta }{1+ϵ\theta }\right)=0,$$

(7)

with corresponding non-dimensional boundary equations

$$\begin{array}{c}\hfill y=-1:f={\beta }_1{e}^{-\alpha \theta }{f}^{\prime },{\theta }^{\prime }=-Bi\theta ,\\ \hfill y=1:f={\beta }_2{e}^{-\alpha \theta }{f}^{\prime },{\theta }^{\prime }=-Bi\theta ,\end{array}$$

(8)

where $\Omega$ is the third grade fluid material parameter, P is the fluid pressure, G is the pressure gradient, $\alpha$ is the variable viscosity parameter, $ϵ$ is the activation energy parameter, $Gr$ is the Grashof number, $Br$ is the Brinkman number, ${\beta }_1$, ${\beta }_2$ are the lower and upper wall slip parameters respectively, $Bi$ is the biot number, $Da$ is the Darcy number, $\delta$ is the porous medium shape parameter, $\lambda$ is the Frank-Kamenetsikii parameter, $\theta$ is the non-dimensional fluid temperature, x and y are the non-dimensional axial and normal coordinates to the inclined channel and f is the non-dimensional fluid velocity.

3. Entropy Generation Analysis

The local volumetric rate of entropy generation for a viscous incompressible third grade fluid in the presence of porous media will be expressed as

$$E_G=\frac{k}{T_0}{\displaystyle \frac{d^{2}T}{d{\overline{y}}^2}}+\frac{\overline{\mu }\left(T\right)}{T_0}{\left({\displaystyle \frac{d\overline{u}}{d\overline{y}}}\right)}^2\left(1+\frac{2{\beta }_3}{\overline{\mu }\left(T\right)}{\left({\displaystyle \frac{d\overline{u}}{d\overline{y}}}\right)}^2\right)+{\displaystyle \frac{\overline{\mu }\left(T\right){\overline{u}}^2}{K{T}_0}},$$

(9)

where the first term is the heat transfer irreversibility, the second term is the entropy generation due to third grade fluid viscous dissipation, and the last term is the irreversibility due to the presence of porous media. Using the non−dimensional parameters and variables (5), we get the equation for the non−dimensional entropy generation rate as

$$N_s={\theta }^{\prime \prime }+Br{e}^{-\alpha \theta }\left[f^{\prime 2}(1+2\Omega e^{\alpha \theta }{f}^{\prime 2})+{\delta }^2{f}^2\right].$$

(10)

The non-dimensional entropy generation number, $N_s=\frac{E_G{a}^{2}E}{kR{T}_0}$, is the ratio of the volumetric entropy generation rate to the characteristic entropy generation rate. Equation (10) can be expressed alternatively as

$$\begin{array}{c}\hfill N_s=N_h+N_f,\\ \hfill N_h={\theta }^{\prime \prime },\\ \hfill N_f=Br{e}^{-\alpha \theta }\left[f^{\prime 2}(1+2\Omega e^{\alpha \theta }{f}^{\prime 2})+{\delta }^2{f}^2\right],\end{array}$$

(11)

where $N_h$ is the entropy generation due to heat transfer irreversibility and $N_f$ is the entropy generation due to the combined fluid viscous dissipation and porosity of the porous medium. The Bejan number $\left(Be\right)$ is used to evaluate the irreversibility distribution and defined as

$$Be=\frac{N_h}{N_h+N_f}.$$

(12)

According to Bejan [19], the Bejan number varies from 0 to 1, where $Be=0$ tells us that the irreversibility due fluid viscous dissipation and porosity dominate, whereas, $Be=\frac{1}{2}$ indicates that the irreversibility due fluid viscous dissipation and porosity is equal to irreversibility due heat transfer in the entropy production, and $Be=1$ reveals that irreversibility due to heat transfer dominates.

4. Method of Solution

The system of nonlinear ordinary differential Equations (6)–(7), with boundary conditions (8) is solved numerically by the spectral quasilinearisation method (SQLM). Firstly, quasilinearisation iteration scheme is applied to linearize the system. Bellman and Kalaba [20] introduced this quasilinearisation method (QLM) which is a generalisation of the Newton-Raphson method for solving nonlinear boundary value problems. The quasilinearisation technique is utilised to linearise the equations before they are solved using the SQLM. We let

$$F=e^{-\alpha \theta }{f}^{\prime \prime }+6\Omega f^{\prime \prime }{f}^{\prime 2}-\alpha e^{-\alpha \theta }{f}^{\prime }{\theta }^{\prime }-\delta e^{-\alpha \theta }f+Gr\theta Sin\left(\omega \right)+G=0,$$

(13)

$$H={\theta }^{\prime \prime }+Br{f}^{\prime 2}(e^{-\alpha \theta }+2\Omega f^{\prime 2})+\delta Br{e}^{-\alpha \theta }{f}^2+\lambda {(1+ϵ\theta )}^{m}exp\left(\frac{\theta }{1+ϵ\theta }\right)=0.$$

(14)

Reducing the non-linear Equations (13) and (14) to linear equations using QLM, we let $f_r$, ${\theta }_r$ be an approximate current solution and $f_{r+1}$, ${\theta }_{r+1}$ be an improved solution of the system of equations. Assuming $|f_{r+1}-f_r|\ll 1$ and $|{\theta }_{r+1}-{\theta }_r|\ll 1$, we linearise the equations by expanding F and H using Taylor series expansion. We have

$$a_{2,r}{f}_{r+1}^{\prime \prime }+a_{1,r}{f}_{r+1}^{\prime }+a_{0,r}f+b_{1,r}{\theta }_{r+1}^{\prime }+b_{0,r}{\theta }_{r+1}=R_{1,r},$$

(15)

$$c_{2,r}{\theta }_{r+1}^{\prime \prime }+c_{1,r}{\theta }_{r+1}^{\prime }+c_{0,r}\theta +d_{1,r}{f}_{r+1}^{\prime }+d_{0,r}{f}_{r+1}=R_{2,r},$$

(16)

with new non−dimensional boundary conditions

$$\begin{array}{c}\hfill y=-1:f_{r+1}={\beta }_1{e}^{-\alpha {\theta }_r}{f}_r^{\prime },{\theta }_{r+1}^{\prime }=-Bi{\theta }_r,\\ \hfill y=1:f_{r+1}={\beta }_2{e}^{-\alpha {\theta }_r}{f}_r^{\prime },{\theta }_{r+1}^{\prime }=-Bi{\theta }_r,\end{array}$$

(17)

where

$$\begin{array}{c}\hfill a_{2,r}={\displaystyle \frac{\partial F}{\partial f_r^{\prime \prime }}}=e^{-\alpha {\theta }_r}+6\Omega f_r^{\prime 2},\\ \hfill a_{1,r}={\displaystyle \frac{\partial F}{\partial f_r^{\prime }}}=-\alpha e^{-\alpha {\theta }_r}{\theta }_r^{\prime }+12\Omega f_r^{\prime \prime }{f}_r^{\prime },\\ \hfill a_{0,r}={\displaystyle \frac{\partial F}{\partial f_r}}=-\delta e^{-\alpha {\theta }_r},\\ \hfill b_{1,r}={\displaystyle \frac{\partial F}{\partial {\theta }_r^{\prime }}}=-\alpha e^{-\alpha {\theta }_r}{f}_r^{\prime },\\ \hfill b_{0,r}={\displaystyle \frac{\partial F}{\partial {\theta }_r}}=-\alpha e^{-\alpha {\theta }_r}{f}_r^{\prime \prime }+{\alpha }^2{e}^{-\alpha {\theta }_r}{f}_r^{\prime }{\theta }_r^{\prime }-\alpha \delta e^{-\alpha {\theta }_r}{f}_r+GrSin\left(\omega \right),\\ \hfill c_{2,r}={\displaystyle \frac{\partial H}{\partial {\theta }_r^{\prime \prime }}}=1,\\ \hfill c_{1,r}={\displaystyle \frac{\partial H}{\partial {\theta }_r^{\prime }}}=0,\\ \hfill c_{0,r}={\displaystyle \frac{\partial H}{\partial {\theta }_r}}=-\alpha Br{e}^{-\alpha {\theta }_r}{f}_r^{\prime 2}-\alpha \delta Br{e}^{-\alpha {\theta }_r}{f}_r^2+\lambda {(1+ϵ\theta )}^{m-2}exp\left(\frac{\theta }{1+ϵ\theta }\right)(mϵ(1+ϵ\theta )+1),\\ \hfill d_{1,r}={\displaystyle \frac{\partial H}{\partial f_r^{\prime }}}=2Br{e}^{-\alpha {\theta }_r}{f}_r^{\prime }+8Br\Omega f_r^{\prime 3},\\ \hfill d_{0,r}={\displaystyle \frac{\partial H}{\partial f_r^{\prime }}}=2\delta Br{e}^{-\alpha {\theta }_r}{f}_r,\\ \hfill R_{1,r}=e^{-\alpha {\theta }_r}{f}_r^{\prime \prime }-\alpha e^{-\alpha {\theta }_r}{\theta }_r^{\prime }{f}_r^{\prime }+6\Omega f_r^{\prime \prime }{f}_r^{\prime 2}-\delta e^{-\alpha {\theta }_r}{f}_r,\\ \hfill R_{2,r}=Br{f}_r^{\prime 2}(e^{-\alpha {\theta }_r}+2\Omega f_r^{\prime 2})+\delta Br{e}^{-\alpha {\theta }_r}{f}_r^2+\lambda {(1+ϵ{\theta }_r)}^{m}exp\left(\frac{\theta }{1+ϵ{\theta }_r}\right).\end{array}$$

(18)

The QLM iteration scheme (15) and (16) is solved using the Chebyshev spectral collocation method. The approximations for the unknown functions are done by using Chebyshev interpolating polynomials in such a way that they are collocated at the Gauss-Lobatto points defined as

$${\tau }_i=cos\left(\frac{\pi i}{N}\right),i=0,1,2,...,N,$$

(19)

where N is the number of collocation points. The functions $f_{r+1}$ and ${\theta }_{r+1}$ at the collocation points are represented by

$$\begin{array}{c}\hfill f_{r+1}\left(\tau \right)=\sum _{k=0}^N{f}_{r+1}\left({\tau }_k\right)T_k\left({\tau }_i\right),\\ \hfill {\theta }_{r+1}\left(\tau \right)=\sum _{k=0}^N{\theta }_{r+1}\left({\tau }_k\right)T_k\left({\tau }_i\right),i=0,1,2,...,N,\end{array}$$

(20)

where $T_k$ is the $k^{th}$ Chebyshev polynomial defined by

$$T_k\left(\tau \right)=cos\left(kco{s}^{-1}\tau \right).$$

(21)

The derivatives of $f_{r+1}$ and ${\theta }_{r+1}$ at the collocation points are represented as

$$\begin{array}{c}\hfill \frac{d^p{f}_{r+1}}{d{y}^p}=\sum _{k=0}^N{\mathbf{D}}_{ki}^p{f}_{r+1}\left({\tau }_k\right),\\ \hfill \frac{d^p{\theta }_{r+1}}{d{y}^p}=\sum _{k=0}^N{\mathbf{D}}_{ki}^p{\theta }_{r+1}\left({\tau }_k\right),\end{array}$$

(22)

where p is the order of differentiation and $\mathbf{D}$ is the Chebyshev spectral differential matrix of order $(N+1)\times (N+1)$. We substitute Equations (19)–(22) into Equations (15)–(16) to obtain

$$\left[a_{2,r}{\mathbf{D}}^2+a_{1,r}\mathbf{D}+a_{0,r}\mathbf{I}\right]{\mathbf{F}}_{\mathbf{r}+\mathbf{1}}+\left[b_{1,r}\mathbf{D}+b_{0,r}\mathbf{I}\right]{\Theta }_{\mathbf{r}+\mathbf{1}}={\mathbf{R}}_{\mathbf{1},\mathbf{r}},$$

(23)

$$\left[d_{1,r}\mathbf{D}+d_{0,r}\mathbf{I}\right]{\mathbf{F}}_{\mathbf{r}+\mathbf{1}}+\left[c_{2,r}{\mathbf{D}}^2+c_{1,r}\mathbf{D}+c_{0,r}\mathbf{I}\right]{\Theta }_{\mathbf{r}+\mathbf{1}}={\mathbf{R}}_{\mathbf{2},\mathbf{r}},$$

(24)

where $\mathbf{I}$ is $(N+1)\times (N+1)$ identity matrix.

Putting Equations (23) and (24) in a matrix form we get

$$\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\left[\begin{array}{c}{\mathbf{F}}_{\mathbf{r}+\mathbf{1}}\\ {\Theta }_{\mathbf{r}+\mathbf{1}}\end{array}\right]=\left[\begin{array}{c}{\mathbf{R}}_{\mathbf{1},\mathbf{r}}\\ {\mathbf{R}}_{\mathbf{2},\mathbf{r}},\end{array}\right]$$

(25)

where

$$\begin{array}{c}\hfill A_{11}=a_{2,r}{\mathbf{D}}^2+a_{1,r}\mathbf{D}+a_{0,r}\mathbf{I},\\ \hfill A_{12}=b_{1,r}\mathbf{D}+b_{0,r}\mathbf{I},\\ \hfill A_{21}=d_{1,r}\mathbf{D}+d_{0,r}\mathbf{I},\\ \hfill A_{22}=c_{2,r}{\mathbf{D}}^2+c_{1,r}\mathbf{D}+c_{0,r}\mathbf{I},\\ \hfill {\mathbf{F}}_{\mathbf{r}+\mathbf{1}}={[f_{r+1}\left({\tau }_0\right),f_{r+1}\left({\tau }_1\right),f_{r+1}\left({\tau }_0\right),...,f_{r+1}\left({\tau }_{i-1}\right),f_{r+1}\left({\tau }_i\right)]}^T,\\ \hfill {\Theta }_{\mathbf{r}+\mathbf{1}}={[{\theta }_{r+1}\left({\tau }_0\right),{\theta }_{r+1}\left({\tau }_1\right),{\theta }_{r+1}\left({\tau }_0\right),...,{\theta }_{r+1}\left({\tau }_{i-1}\right),{\theta }_{r+1}\left({\tau }_i\right)]}^T.\end{array}$$

(26)

The approximate solutions for $\mathbf{F}$ and $\Theta$ are obtained by solving the matrix system (25),

$$\left[\begin{array}{c}{\mathbf{F}}_{\mathbf{r}+\mathbf{1}}\\ {\Theta }_{\mathbf{r}+\mathbf{1}}\end{array}\right]={\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]}^{-1}\left[\begin{array}{c}{\mathbf{R}}_{\mathbf{1},\mathbf{r}}\\ {\mathbf{R}}_{\mathbf{2},\mathbf{r}}\end{array}\right]$$

(27)

and we use a suitable initial approximation that satisfies the boundary conditions of Equations (6) and (7) as

$$f_0\left(y\right)=-\frac{1}{2}G{y}^2-\frac{G({\beta }_1+{\beta }_2)y}{{\beta }_1-{\beta }_2+2}-\frac{1}{2}\frac{G(4{\beta }_1{\beta }_2-3{\beta }_1+3{\beta }_2-2)}{{\beta }_1-{\beta }_2+2},$$

(28)

$${\theta }_0\left(y\right)=0.$$

(29)

5. Results and Discussion

Combined effects of temperature dependent viscosity, porous media permeability, Navier slip and convective boundary conditions on the flow of a reactive third grade fluid in an inclined channel have been studied. In the following discussion, we carry out a detailed thermodynamic analysis of the dependence of the fluid velocity and temperature profiles, the entropy generation and the irreversibility ratio on the various thermophysical parameters embedded in the flow system. We present the SQLM graphical solutions and provide a qualitative description of the simulated variations. Unless otherwise stated, the following parameter values are employed: $\alpha =0.1$, $\Omega =0.1$, $\delta =0.1$, $Gr=0.8$, $\omega =Pi/4$, $G=1$, $Br=0.5$, $\lambda =0.1$, $ϵ=0.1$, ${\beta }_1=0.1$, ${\beta }_2=0.1$, $Bi=5.0$, and $m=0.5$. Where a parameter is not varied, these will be the default values used.

5.1. Residual Analysis

The accuracy and convergence of the SQLM were validated by performing a residual error analysis. The residual errors were calculated to gain understanding of the accuracy, as defined by

$$\begin{array}{c}\hfill {\left|\right|Res\left(f\right)\left|\right|}_{\infty }=\left|\right|e^{-\alpha \theta }{f}^{\prime \prime }+6\Omega f^{\prime \prime }{f}^{\prime 2}-\alpha e^{-\alpha \theta }{f}^{\prime }{\theta }^{\prime }-\delta e^{-\alpha \theta }f+Gr\theta Sin\left(\omega \right)+G\left|\right|,\\ \hfill {\left|\right|Res\left(\theta \right)\left|\right|}_{\infty }=\left|\right|{\theta }^{\prime \prime }+Br{f}^{\prime 2}(e^{-\alpha \theta }+2\Omega f^{\prime 2})+\delta Br{e}^{-\alpha \theta }{f}^2+\lambda {(1+ϵ\theta )}^{m}exp\left(\frac{\theta }{1+ϵ\theta }\right)\left|\right|.\end{array}$$

The error infinity norms, which are used to confirm the convergence of the SQLM approximations are defined by

$$\begin{array}{c}\hfill \mathrm{Error}\left(\mathrm{f}\right)=\underset{0\le i\le N}{max}{\Vert {\mathbf{F}}_{\mathbf{r}+\mathbf{1},\mathbf{i}}-{\mathbf{F}}_{\mathbf{r},\mathbf{i}}\Vert }_{\infty },\\ \hfill Error\left(\theta \right)=\underset{0\le i\le N}{max}{\Vert {\Theta }_{\mathbf{r}+\mathbf{1},\mathbf{i}}-{\Theta }_{\mathbf{r},\mathbf{i}}\Vert }_{\infty }.\end{array}$$

We note that the accuracy of order ${10}^{-5}$ for residual errors were attained after two iterations for $f\left(y\right)$ and $\theta \left(y\right)$ as shown in Table 1, and this shows that the method is highly accurate with a good convergence rate. To have a clearer picture on the convergence rates, we plot the residual errors against the number of iterations in Figure 2. Figure 2 shows that an increase in the number of iteration results in a decrease in the error infinity norm and the method converge after two iterations.

5.2. Velocity Profile

The response of the fluid velocity profile to the various flow parameters is displayed in Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. The figures reveal parabolic velocity profiles where the maximum velocity is recorded at the centre of the channel. Figure 3 reveals that from the centre of the channel towards the lower wall, the velocity profile decreases with increasing values of the variable viscosity parameter $\alpha$. On the other hand, towards the upper channel wall, the figure shows the profile increasing with the variable viscosity parameter. Bearing in mind that increasing values of the variable viscosity parameter means that the fluid viscosity is reduced, this observation is to be expected since the lower wall is subjected to convective heat exchange with the ambient and on the upper wall there is a constant heat flux. As the fluid emits heat into the ambient the fluid viscosity increases resulting in velocity retardation. In Figure 4 the Brinkmann number $Br$ is observed to affect the fluid velocity in a way similar to the variable viscosity parameter, although its effect is less pronounced as compared to the effect of the variable viscosity parameter. Since the Brinkmann number is a measure of heat flux from the channel plates to the viscous fluid, the observed similarity is to be expected.

In Figure 5 we see that the activation energy parameter $ϵ$ has no significant effect on the fluid velocity. This is explained by the fact that the activation energy parameter only enters the velocity equation implicitly through the temperature/viscosity coupling. In this way, its effects on the velocity can at best be marginal. An increase in the Grashof number $Gr$, Figure 6, is observed to increase the fluid velocity. Variations in the Grashof number explain the buoyancy effects on the flow. Higher Grashof numbers point to increased buoyancy source terms that result in the corresponding increase in the fluid velocity as illustrated in Figure 6. The variation of the fluid velocity in response to increasing values of the Frank-Kamenetsikii parameter (chemical reaction parameter $\lambda$) is displayed in Figure 7. An increase in the chemical reaction parameter means that the rate of chemical reaction increases, and so are the heating source terms. Figure 7 shows that towards the lower channel wall, the reaction parameter retards the velocity, whereas towards the upper wall the velocity is enhanced. At the lower wall, as explained earlier, convective heat emitted out into the ambient means that the reaction rate is suppressed and so the fluid particles have less energy to move. At the upper wall, a constant heat flux into the fluid boosts the rate of chemical reaction that results in increased excitation of the fluid particles.

Figure 8 and Figure 9 show the flow damping effects of the third grade material parameter and the porous medium shape parameter $\delta$, respectively. An increase in the third grade material parameter $\Omega$ results in corresponding increases in the non-Newtonian properties of the fluid, like the visco-elastic effects and others, that induce increased resistance to flow. On the other hand, increasing the porous medium shape parameter means that the pore spaces in the porous matrix are reduced, resulting in reduced porosity and hence damping the flow.

5.3. Fluid Temperature Distribution

Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16 display variations of the fluid temperature profile in response to various flow governing parameters. A common feature that is revealed by all these graphs is the dipping of the fluid temperature to very low levels as we move closer to and at the lower channel wall. This is attributed to the high Biot number ($Bi=5.0$) at the lower channel wall that was set as the default value in all the computations. In general, high values of the Biot number render increased convective heat loss from the wall to the ambient, leading to a decrease in the wall temperature. Higher degrees of convective cooling at the lower wall induces a significant temperature drop in the bulk of the fluid closer to the wall.

Except for Figure 12, where the activation energy parameter shows only marginal effects on the fluid temperature, and Figure 15, where the third grade material parameter retards the fluid temperature, all the other figures show the fluid temperature increasing with an increase in the values of the variable viscosity parameter, the Brinkmann number, the Grashof number, the chemical reaction parameter, and the porous medium shape parameter. The Brinkmann number and the chemical reaction parameter are observed to cause the most significant increase in fluid temperature. As alluded to earlier in the discussion, increasing values of the variable viscosity parameter means that the fluid viscosity is reduced, and this diminishes the fluid’s resistance to flow. The resultant increased fluid velocity increases the viscous heating source terms in the temperature equation, leading to the elevation of fluid temperature as depicted in Figure 10. In Figure 11, the Brinkmann number, a measure of heat flux from the channel wall into the fluid, is observed to increase the fluid temperature significantly as pointed to earlier. In Figure 13, increasing buoyancy forces lead to an increase in fluid temperature. The heat flux from the upper channel wall into the fluid leads to a volumetric thermal expansion of the fluid which induces some overturning instability in the fluid that causes the fluid particles to move faster. Due to coupling, the increased velocity increases the heating source terms in the temperature equation, resulting in enhanced fluid temperature. Figure 14 shows the fluid temperature increasing significantly with an increase in the values of the chemical reaction parameter. This phenomenon has been explained earlier in Section 5.2

In Section 5.2, we explained the reduction of fluid velocity due to increasing non-Newtonian character as measured by the third grade material parameter. The reduced velocity in turn decreases the viscous heating source terms in the temperature equation which leads to a drop in fluid temperature as shown in Figure 15. In Figure 16, except very close to the upper wall, elsewhere in the channel an increase in the porous medium shape parameter values cause a slight increase in the fluid temperature. This can be attributed to the frictional forces caused by the fluid particles as they force their way into reduced pore spaces. Close to the upper wall, reduced porosity is seen to decrease the fluid temperature. This, again, can be explained by the coupling effect via damped velocity.

5.4. Entropy Generation

Analysis of the entropy generated in the flow is carried out with the aid of Figure 17, Figure 18, Figure 19, Figure 20, Figure 21, Figure 22, Figure 23. The figures depict parabolic entropy profiles with minimum entropy rate $Ns=0$ at the channel core region. Except for Figure 20, all the other figures show that the entropy generation at the core region of the channel remains unaffected by the parameter variations. Figure 17, Figure 18, Figure 19 and Figure 20 show the entropy generation rate increasing with increasing values of the Brinkmann number, the activation energy parameter, the Grashof number and the chemical reaction parameter, respectively. The most significant effects are in Figs Figure 17 and Figure 20, due to the Brinkmann number and the chemical reaction parameter. In Figure 20, the rate of entropy generation is increased throughout the entire channel including the core region. In Figure 21 and Figure 22 the third grade material parameter and the porous medium shape parameter have a retardation effect to the entropy generation rate. In the lower half of the channel, Figure 23, an increase in the variable viscosity parameter enhances the entropy generation rate, while in the upper half of the channel the variable viscosity parameter has the opposite effect. Comparing Figure 23 and Figure 3, we observe that in the region where the variable viscosity parameter decreases the fluid velocity, it enhances the rate of entropy generation. In the region where the variable viscosity parameter increases the fluid velocity, it retards the entropy generation rate.

5.5. Bejan Number

Figure 24, Figure 25, Figure 26, Figure 27, Figure 28, Figure 29 and Figure 30 demonstrate the irreversibility ratio trends in response to variation of the flow governing parameters. A striking feature displayed by all these 7 graphs is that the irreversibility ratio is maximum at the core region of the channel. This shows that, in this flow, irreversibility due to heat transfer dominates at the channel core region. Elsewhere within the channel, the graphs reveal the combined irreversibility due to viscous dissipation and porous medium porosity dominating the heat transfer irreversibility at the core of the lower half of the channel (at $y=-0.5$) and at the core of the upper half (at $y=0.5$). Figure 24 shows that as the variable viscosity parameter is increased, from $y=-1$ to $y=-0.5$ the Bejan number increases and in the next quarter of the channel the opposite effect happens. From the centre of the channel to $y=0.5$ the Bejan number increases and in the last quarter of the channel it is retarded. An interesting observation in Figure 25 is that the irreversibility ratio is constant at Be=1 throughout the channel when $Br=0$. The same observation is made in Figure 26 when $ϵ=0$. As the Brinkmann number increases, the Bejan number increases in the lower quarter of the channel in Figure 25. In the same graph, from $y=-0.5$ to $y=0.5$, the Bejan number decreases with the Brinkmann number, and in the upper quarter it is enhanced. In Figure 26, increasing values of the activation energy parameter diminishes the Bejan number in both the lower half and the upper half of the channel.

The effects of the Grashof number on the irreversibility ratio, Figure 27, mirrors those of the Brinkmann number, albeit at a less pronounced intensity. In Figure 28, an increase in the chemical reaction parameter reveals effects that are exactly opposite of the effects of the activation energy parameter described in Figure 26. In Figure 29, an increase in the non-Newtonian properties of the fluid decreases the Bejan number in the lower half of the channel, whereas in the upper half of the channel it is enhanced. Figure 30 shows that increasing the porous medium shape parameter retards the Bejan number at the core region of the channel, enhances it at the centre of both halves of the channel and retards it again as we approach the channel walls.

6. Conclusions

In this article, entropy generation analysis of a steady flow of a reactive variable viscosity third grade fluid through a porous saturated medium with Navier slip and convective boundary conditions was successfully carried out with the aid of SQLM. The numerical method was found to be highly accurate and convergent. Graphical analysis of the response of the velocity and temperature distributions, entropy generation rate and irreversibility ratio to the various parameters embedded in the flow system revealed that:

• The temperature dependent viscosity, the Brinkmann number and the chemical reaction parameter have a retardation effect on the fluid velocity in the lower half of the channel that is subjected to convective heat exchange with the ambient, while the opposite effect is observed in the upper half of the channel that is subjected to a constant heat flux.
• The non-Newtonian properties of the fluid and the porous medium shape factor parameter were observed to have a damping effect on the fluid velocity, while the buoyancy forces were observed to have the opposite effect.
• The fluid temperature profile was observed to increase with an increase in the Brinkmann number, the buoyancy force, and the chemical reaction parameter, whereas the variable viscosity parameter and the third grade material parameter had a retardation effect.
• A parabolic entropy generation rate profile was observed with minimum entropy $N_s=0$ at the channel core region.
• Parameters either increase or decrease the entropy generation rate elsewhere in the channel but at the core region of the channel it remained unaffected by all the parameters except the chemical reaction parameter.
• The irreversibility ratio was observed to be maximum at the channel core region indicating dominance of heat transfer irreversibility over the combined viscous dissipation and porosity irreversibility.
• Elsewhere in the channel parameters revealed varied effects on the Bejan number but in general irreversibility due to viscous dissipation and porosity was found to dominate heat transfer irreversibility at the core region of either half of the channel.

These observed interesting manifestations will no doubt inform optimal system regulation criteria for entropy generation minimization in thermomechanical systems to achieve energy utilization efficiency. There is growing demand for economical and maintainable systems across the globe. Most industrial and engineering flow processes and thermal systems are unable to work at optimal level due to entropy generation. Since entropy generation is a measure of the destruction of accessible work of the systems, the determination of entropy generation is thus extremely vital for upgrading system performance. Our study is a theoretical model through which the factors (or parameters) that contribute to entropy generation can be identified so that their effects can be minimized through intelligent regulation so that the flow system efficiency can be maximised. Moreover, with the global technological advancements on non-Newtonian fluids, the theoretical results will provide guidance on improving, conserving and upgrading several designs in industrial thermo-fluid systems.