Nonlinear Mixed Convection in a Reactive Third-Grade Fluid Flow with Convective Wall Cooling and Variable Properties

Abstract

Energy management and heat control whenever a reactive viscous fluid is the working medium has been one of the greatest challenges encountered by many in the field of chemical and industrial engineering. A mathematical approach to thedetermination of critical points beyond which the working environment becomes hazardous is presented in the present investigation together with the entropy generation analysis that guarantees the efficient management of expensive energy resources. In this regard, the nonlinear mixed convective flow behavior of a combustible third-grade fluid through a vertical channel with wall cooling by convection is investigated. The mathematical formulation captures the nonlinearities arising from second-order Boussinesq approximation and exponential dependence of internal heat generation, viscosity, and thermal conductivity on temperature. The resulting nonlinear boundary value problems were solved based on the spectral Chebyshev collocation method (SCCM) and validated with the shooting-Runge–Kutta method (RK4). The nonlinear effects on the flow velocity, temperature distribution, entropy generation, and Bejan heat irreversibility ratio are significant. Further analyses include the thermal stability of the fluid. Findings from the study revealed that flow, temperature, and entropy generation are enhanced byincreasing values of the Grashof number, the quadratic component of buoyancy, and the Frank-Kameneskii parameter, but are reducedbyincreasing the third-grade material parameter. Moreover, it was shown that increasing values of the third-grade parameter encourages the thermal stability of the flow, while increasing values of the linear and nonlinear buoyancy parameter destabilizes the flow. The present result is applicable to thick combustible polymers with increased molecular weight.

1. Introduction

A renewed interest in energy utilization and management has led to a significant increase in studies centered on a broad range of reactive non-Newtonian fluids in the last few decades. A good number of researchers gave mathematical explanations from a fluid dynamics viewpoint on how best to maximize exergy when working with reactive fluids. For instance, Salawu et al. [1] reported the transient flow of a reactive fourth-order fluid through a horizontal channel

Permeated by a magnetic field, Cuiet al.in [2] investigated the two-dimensional reacting flow of hydromagnetic nano-Oldroyd-B fluid driven by mixed convection. In [3], Salawu and his cohorts considered the slip flow of combustible MHD couple stress fluids with variable properties. Sadiq and Hayat [4] presented a model for optimizing entropy generation in Reiner–Rivlin flow over a stretched rotating disc. Okoya [5] reported a numerical investigation of a time-independent variable viscous third-grade fluid flow through a cylinder. Adesanya et al. [6] analyzed the flow of third-grade fluid between two solid boundaries under Arrhenius kinetics. The study in [7] focused on the reactive flow under the Erying–Powell constitutive model. More studies on reactive flows of third-grade fluid can be found in [8,9,10,11,12,13,14] and the literature cited therein.

The quantity of energy absorbed by the flow is defined by the degree of friction caused by viscosity, which plays a crucial role in the creation of fluid measures for flow measurements. The changing viscosity in Couette and Poiseuille flows of a tangential hyperbolic fluid via an inclined channel with porous media was explored by Zehra et al. [15]. The fluid’s viscosity has been assumed in this research to be inversely proportional to a linear function of pressure. It should go without saying that the viscosity and thermal conductivity of any fluid may be alteredas the temperature rises. This phenomenon therefore plays a crucial role in the rate of heat transport near the surface, much as other thermophysical features. It is acknowledged that certain characteristics, particularly fluid viscosity, may vary with temperature. This change in viscosity and thermal conductivity must be taken into consideration in order to correctly forecast the flow and heat transfer rates. Jeffrey fluid flow, heat transfer, and mass transfer via a permeable wave-like channel with variable viscosity and thermal conductivity were studied by Manjunatha et al. in their study [16]. Qasim et al.’s [17] investigation on the impact of dissipation of a fluid flow through a tiny needle takes these fluid parameters into account. Recently, Saraswathy et al. [18] investigated these phenomena in order to use a numerical approach to investigate the solution for an asymmetric flow and heat transfer of a rotating fluid with the influence of Arrhenius energy.

Additionally, while the first law of thermodynamics has been used by certain writers, its efficacy has been shown to be inferior to that of the second law. Thermal efficiency decreases and system entropy increases as a result of energy losses in engineering and thermal equipment. It is well known that the development of entropy during any thermal process determines the degree of irreversibility. The optimization procedures have recently paid a lot of attention to this examination. Many thermal systems may include irreversible processes. The second law of thermodynamics has been extensively employed in the literature to enhance these types of irreversibility. In this light, Khan et al. [19] analyzed the formation of entropy on an Erying–Powell liquid flowing via a permeable channel using a homotopy analysis technique. Singh et al. [20] investigated the entropy generation impact in a micropolar fluid flowing through an aligned tube. Here, it is claimed that the fluid’s viscosity and thermal conductivity both change with temperature. The conclusion of their analysis shows that temperature-dependent thermal and viscous conductivity factors enhance the entropy production profile. Similar studies on entropy production rates and irreversibility in convection-based heat transport may be found in references [21,22,23].

In modeling the flow driving mechanism, the nonlinear temperature-density relationship has been established over the years to give a better approximation than the linearized form of heat transfer in a wide range of buoyancy-driven flows due to the considerable difference in the fluid temperature and the ambient. Given these numerous applications, Adesanyaet al. [24] focused on the entropy generation, flow, and nonlinear heat transfer to a couple stress fluid flow through aporous medium. In the Xiaet al. [25] report, a mathematical explanation for the nonlinear bio-convective flow of hydromagnetic fluid containing microorganisms with Hall effect is discussed. Ibrahim and Gadisa [26] presented a finite element residual approximation for Oldroyd-B flow with internal heat generation. Patilet al. [27] constructed the numerical approximation to nonlinear convective flow problems with Brownian and thermophoresis. Waqaset al. [28] reported the steady developing flow of Williamson fluid with variable diffusion and thermal conductivity. A semi-analytical approach is employed by Yusuf et al. [29] to examine the magnetohydrodynamic (MHD) nonlinear convective flow of a reactive non-Newtonian fluid model with convective heating. Other vital studies on nonlinear convective flows include refs. [30,31,32] and the cited references.

Despite the huge body of literature surrounding heat irreversibility in a combustible non-Newtonian fluid flow, the nonlineardependence of some fluid parameters on temperatureremains poorly understood with respect to thermal stability for the determination of safe working conditions and energy management or efficiency. The proposed nonlinear coupled model will ensure that the thermal system is not underdetermined when used to predict buoyancy-induced flow and thermal behavior. The nonlinearity of the model suggests the non-existence of a close-form solution; therefore, the present work focuses on a fast-converging numerical approach totackling the coupled nonlinear convective flow and heat transfer problem. To the best of our knowledge, the problem presented here has not been undertaken byany researcher, and that explains the novelty of the work.

Including the nonlinear buoyancy approximation, nonlinear internal heat, and exponential dependence of fluid viscosity and thermal conductivity is expected to produce a better approximation of reactive flow and heat transfer problems, especially during a number of engineering processes, such as advanced powder synthesis by self-propagating high-temperature, heat exchangers, loop/fluidized bed reactors, treatment of reactive hydrocarbons, combustion chambers and many more engineering applications. In the next section, the mathematical analysis of the problem is presented.

2. Mathematical Formulation

In the model formulation, we consider a fully-developed flow of a reactive third-grade fluid flow through a vertical channel with Newtonian cooling at the walls. The buoyancy effect is felt upstream in the vertical x-axis while the y-axis is perpendicular to it, as shown in Figure 1. All the fluid properties like viscosity, thermal conductivity, and internal heat generation are assumed to vary with temperature in an exponential manner, except density which has quadratic dependence.

Under the flow assumptions, the vector form for the momentum equation becomes:

$$\rho \frac{DV}{Dt}=\rho g+div\tau$$

(1)

where the operator,$\frac{D(.)}{Dt}=\left(\frac{\partial }{\partial t}+V.\nabla \right)(.)$ is the material derivate acting on the velocity vector $V$, and $\rho ,g$-density and gravitational force, respectively. The last term $\tau$ is the Cauchy stress tensor defined as [14,33]

$$\tau =-PI+\mu S_1+{\alpha }_1{S}_2+{\alpha }_2{S}_1^2+{\beta }_1{S}_3+{\beta }_2\left(S_1{S}_2+S_2{S}_1\right)+{\beta }_3\left(tr{S}_1^2\right)S_1$$

(2)

In (2), ${\alpha }_n, \left(n=1,2,3\right),{\beta }_n, \left(n=1,2,3\right)$ are material effects, $\left(P,I,\mu \right)$ are the pressure, tensor identity and fluid viscosity, respectively while $S_n, \left(n=1,2,3\right)$ represents the Rilvin–Ericksen tensors and are defined as

$$\begin{array}{l}{S}_1=\nabla V+{\left(\nabla V\right)}^T\\ S_n=\frac{D{S}_{n-1}}{Dt}+S_{n-1}\nabla V+{\left(\nabla V\right)}^T{S}_{n-1}, n>1\end{array}$$

(3)

The Clausius–Duhem inequality condition $\mu \ge 0, {\alpha }_1\ge 0 ,\left|{\alpha }_1+{\alpha }_2\right|\le \sqrt{24\mu {\beta }_3}$ ensures the thermodynamic process is consistent when ${\beta }_1={\beta }_2=0$ implies that ${\beta }_3\ge 0$.

Therefore, (2) reduces to

$$\tau =-PI+\left(\mu +{\beta }_3\left(tr{S}_1^2\right)\right)S_1+{\alpha }_1{S}_2+{\alpha }_2{S}_1^2$$

(4)

From (4), it is easy to note that the effective viscosity, which is assumed to be a nonlinear function of temperature, becomes

$${\mu }_{eff}=\mu \left(T\right)+{\beta }_3\left(tr{S}_1^2\right)$$

(5)

As a consequence of (5), the momentum equation with the nonlinear temperature-dependent fluid density becomes

$$0=-\frac{dP}{dx}+\frac{d}{dy{}^{\prime }}\left(\mu \left(T\right)\frac{du{}^{\prime }}{dy{}^{\prime }}+2{\beta }_3{\left(\frac{du{}^{\prime }}{dy{}^{\prime }}\right)}^3\right) +\rho g{\varphi }_0\left(T-T_0\right)+\rho g{\varphi }_1{\left(T-T_0\right)}^2$$

(6)

While the balanced energy equation with variability in thermal conductivity becomes

$$0=\frac{d}{dy{}^{\prime }}\left(k\left(T\right)\frac{dT}{dy{}^{\prime }}\right)+Q{C}_{0}A{\left(\frac{hT}{\upsilon l}\right)}^m{e}^{-\frac{E}{RT}} +\left(\mu \left(T\right)+2{\beta }_3{\left(\frac{du{}^{\prime }}{dy{}^{\prime }}\right)}^2\right){\left(\frac{du{}^{\prime }}{dy{}^{\prime }}\right)}^2$$

(7)

The fact that the heat transfer to the viscous non-Newtonian fluid is irreversible suggests that the entropy generation profile is a nonlinear combination of heat transfer and fluid friction, that is,

$$E_G=\frac{k\left(T\right)}{T_0^2}{\left(\frac{dT{}^{\prime }}{dy{}^{\prime }}\right)}^2+\frac{1}{T_0}{\left(\frac{du{}^{\prime }}{dy{}^{\prime }}\right)}^2\left(\mu \left(T\right)+2{\beta }_3{\left(\frac{du{}^{\prime }}{dy{}^{\prime }}\right)}^2\right)$$

(8)

Togetherwith the essential and temperature dependent natural boundary conditions

$$\begin{array}{l}u{}^{\prime }=0, -k\left(T\right)\frac{dT}{dy{}^{\prime }}=h_1\left(T_0-T\right) y{}^{\prime }=0\\ u{}^{\prime }=0, -k\left(T\right)\frac{dT}{dy{}^{\prime }}=h_2\left(T-T_0\right) y{}^{\prime }=a\end{array}$$

(9)

Of engineering interest is the wall skin friction and Nusselt number for heat transfer at the wall

$$C_f={{\mu }_0{e}^{-\overline{\alpha }\left(T-T_0\right)}\frac{du{}^{\prime }}{dy{}^{\prime }}+2{\beta }_3 {\left(\frac{du{}^{\prime }}{dy{}^{\prime }}\right)}^3|}_{y{}^{\prime }=0}, q_w= {k_0{e}^{-\overline{\eta }\left(T-T_0\right)}\frac{dT}{dy{}^{\prime }}|}_{y{}^{\prime }=0}$$

(10)

The Reynold’s exponential dependence of viscosity and thermal conductivity on temperature is assumed to take the form

$$\mu \left(T\right)={\mu }_0{e}^{-\overline{\alpha }\left(T-T_0\right)}, k\left(T\right)=k_0{e}^{-\overline{\eta }\left(T-T_0\right)}$$

(11)

With the introduction of the following variables and parameters

$$\begin{array}{l}y=\frac{y{}^{\prime }}{a}, u=\frac{u{}^{\prime }}{UG}, \theta =\frac{E\left(T-T_0\right)}{R{T}_0^2},\alpha =\frac{\overline{\alpha }R{T}_0^2}{E},\gamma =\frac{{\beta }_3{U}^2{G}^2}{{\mu }_0{a}^2},\delta ={\left(\frac{\upsilon \iota }{\hslash T_0}\right)}^m\frac{{\mu }_0{U}^2{G}^2}{QA{C}_0{a}^2}{e}^{\frac{E}{R{T}_0}}\\ G=-\frac{a^2}{{\mu }_{0}UG}\frac{dP}{dx}=1,Gr=\frac{\rho g{\varphi }_0{a}^{2}R{T}_0^2}{{\mu }_{0}EUG},\sigma =\frac{{\varphi }_{1}R{T}_0^2}{{\varphi }_{0}E},\lambda ={\left(\frac{\hslash T_0}{\upsilon \iota }\right)}^m\frac{QEA{a}^2{C}_0}{k_{0}R{T}_0^2}{e}^{-\frac{E}{R{T}_0}},B{i}_{1,2}=\frac{h_{1,2}a}{k_0}\end{array}\}$$

(12)

We obtainthe following dimensionless, nonlinear, coupled, system of boundary value problems

$$\begin{array}{l}0=1+\frac{d}{dy}\left(e^{-\alpha \theta }\frac{du}{dy}\right)+6\gamma {\left(\frac{du}{dy}\right)}^2\frac{d^{2}u}{d{y}^2}+Gr\left(\theta +\sigma {\theta }^2\right), u\left(0\right)=0=u\left(1\right);\\ \\ 0=\frac{d}{dy}\left(e^{-\eta \theta }\frac{d\theta }{dy}\right)+\lambda \left({\left(1+\epsilon \theta \right)}^m{e}^{\frac{\theta }{1+\epsilon \theta }}+\delta {\left(\frac{du}{dy}\right)}^2\left(e^{-\alpha \theta }+2\gamma {\left(\frac{du}{dy}\right)}^2\right)\right), \\ {\left(e^{-\eta \theta }\frac{d\theta }{dy}-B{i}_1\theta \right)|}_{y=0} =0 ={\left(e^{-\eta \theta }\frac{d\theta }{dy}+B{i}_2\theta \right)|}_{y=1}. \end{array}\}$$

(13)

Please note that the special case of the coupled Equation (13) is the non-reactive Newtonian flow when $\alpha =\gamma =Gr=\eta =0$, that was reported numerically in [33] by using shooting-Runge–Kutta scheme within Maple environment. In its dimensionless form, the expression for the positive definite entropy generation becomes

$$N_S= e^{-\eta \theta }\underset{N_h}{\underbrace{{\left(\frac{d\theta }{dy}\right)}^2}}+\underset{N_f}{\underbrace{\frac{\lambda \delta }{\epsilon }\left({\left(\frac{du}{dy}\right)}^2\left(e^{-\alpha \theta }+2\gamma {\left(\frac{du}{dy}\right)}^2\right)\right)}}$$

(14)

Component-wise, the Bejan irreversibility ratio becomes

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

(15)

While the expressions for the dimensionless skin friction and Nusselt number becomes

$$S_f=\frac{\rho a^2{C}_f}{{\mu }_0^2}={e^{-\alpha \theta }\frac{du}{dy}+2\gamma {\left(\frac{du}{dy}\right)}^3|}_{y=0}, Nu=\frac{aE{q}_w}{k_{0}R{T}_0^2}= {e^{-\eta \theta }\frac{d\theta }{dy}|}_{y=0}$$

(16)

3. Solution by Weighted Residual Method

A collocation method based on Chebyshev polynomials is considered here to obtain numerical solutions of a differential equationin [0, L]. This technique approximatesa solution to a differential equation

$$\begin{array}{l}f(\varphi (x))+g(x)=0, \mathrm{in} \mathrm{domain}, D=[0,L]\hfill \end{array}$$

(17)

The approximate solution needs to satisfy the boundary conditions and it is given as

$$\begin{array}{l}\psi (x)={\displaystyle \sum _{i=0}^N}{k}_i{\stackrel{⃛}{T}}_i\left(\frac{2x}{L}-1\right)\hfill \end{array}$$

(18)

where $\varphi (x)$ is the unknown dependent function, f(x) is the source term, $k_i$ are undetermined coefficients, ${\stackrel{⃛}{T}}_i\left(\frac{2x}{L}-1\right)$ are shifted Chebyshev base functions from [−1,1] to [0, L]. To acquire the values of $k_i$, Equation (18) is inserted into Equation (17) to yield residual error $R_{\psi }(x,k_i)$. $R_{\psi }(x,k_i)$ is forced to be closer to zero using the collocation method as follows;

$$\begin{array}{l}\mathrm{for} \delta \left(x-x_j\right)=\{\begin{array}{cc}1,& x=x_j\\ 0,& \mathrm{otherwise},\end{array}\hfill \end{array}$$

$$\begin{array}{l}{\displaystyle {\int }_0^1}{R}_{\psi }(x,k_i)\delta \left(x-x_j\right)dx=R_{\psi }(x,k_i)=0, \mathrm{for} j=1, 2,\dots N-1\hfill \end{array}$$

(19)

where $x_j=\frac{1}{2}\left(1-\cos \left(\frac{j\pi }{N}\right)\right)$ are shifted Gauss–Lobatto points. Therefore, Equation (14) as well as the equations generated from the boundary conditions constitute a system of algebraic equations that must be solved in order to obtainthe constant coeficient values, $k_i$. In order to obtain the solutions of Equation (13), solutions are assumed for $u(y)$ and $\theta (y)$, respectively, in the form:

$$\begin{array}{l}u\left(y\right)={\displaystyle \sum _{i=0}^N}{a}_i{\stackrel{⃛}{T}}_i\left(2y-1\right), \theta \left(y\right)={\displaystyle \sum _{i=0}^N}{b}_i{\stackrel{⃛}{T}}_i\left(2y-1\right)\hfill \end{array}$$

(20)

where $a_i$ and $b_i$ are constants to be determined, and ${\stackrel{⃛}{T}}_i\left(2y-1\right)$ are the shifted Chebyshev base functions. To obtain the values of $a_i$ and $b_i$, to satisfy the boundary conditions, Equation (20) is substituted into the boundary conditions in Equation (13) and thus have

$${{\displaystyle \sum _{i=0}^N}\stackrel{⃛}{T}{a}_i(2y-1)|}_{y=0}=0={{\displaystyle \sum _{i=0}^N}\stackrel{⃛}{T}{a}_i(2y-1)|}_{y=1},$$

(21)

$$\begin{array}{l}{e^{-\eta \left({\displaystyle \sum _{i=0}^N}\stackrel{⃛}{T}{b}_i(2y-1)\right)}{\displaystyle \sum _{i=0}^N}\frac{d{\stackrel{⃛}{T}}_i}{dy}{b}_i(2y-1)-B{i}_1{\displaystyle \sum _{i=0}^N}\stackrel{⃛}{T}{b}_i(2y-1)|}_{y=0}=0={e^{-\eta \left({\displaystyle \sum _{i=0}^N}\stackrel{⃛}{T}{b}_i(2y-1)\right)}{\displaystyle \sum _{i=0}^N}\frac{d{\stackrel{⃛}{T}}_i}{dy}{b}_i(2y-1)-B{i}_2{\displaystyle \sum _{i=0}^N}\stackrel{⃛}{T}{b}_i(2y-1)|}_{y=1}\hfill \end{array}$$

(22)

Equation (15) is also substituted into Equations (9) and (10), resulting in $re{s}_u(y,a_i)$ and $re{s}_{\theta }(y,a_i,b_i)$, respectively, as residues. At this juncture, the residues are forced to be zero by applying collocation method as follows

$$\begin{array}{l}\begin{array}{l}{\displaystyle {\int }_0^1}re{s}_u\delta \left(y-y_j\right)dy=re{s}_u\left(y_j,a_i\right)=0, \mathrm{for} j=1,2,\dots N-1\hfill \end{array}\\ \begin{array}{l}{\displaystyle {\int }_0^1}re{s}_{\theta }\delta \left(y-y_j\right)dy=re{s}_{\theta }\left(y_j,a_i,b_i\right)=0, \mathrm{for} j=1,2,\dots N-1\hfill \end{array}\end{array}$$

(23)

where $y_j=\frac{1}{2}\left(1-\cos \left(\frac{j\pi }{N}\right)\right)$. Thus, Equations (16)–(18) form 2N + 2 system of algebraic equations containing 2N + 2 undetermined coefficients $(a_i,b_i)$. This procdeure is coded in a mathematicalsymbolic software to solve the generated equations using the Newton technique.

4. Solution by Shooting-Runge–Kutta Method

To achieve a numerical solution by the Runge–Kutta scheme, we first apply the shooting technique to the coupled problem by converting the second-order differential equations to a set of first-order ordinary Lipshitz continuous differential equations.

$$z_1=u, z_2=u{}^{\prime }, z_3=\theta , z_4=\theta {}^{\prime }$$

(24)

With (8), the coupled BVP becomes a coupled system of first order differential equations

$$\begin{array}{l}{\left(\begin{array}{l}{z}_1\\ \\ z_2\\ \\ z_3\\ \\ \\ z_4\end{array}\right)}^{{}^{\prime }}=\left(\begin{array}{l}{z}_2\\ \frac{\alpha e^{-\alpha z_3}-Gr\left(z_3+\sigma z_3^2\right)-1}{e^{-\alpha z_3}+6\gamma z_2^2}\\ z_4\\ e^{\eta z_3}\left(\eta e^{-\eta z_3}{z}_4^2-\lambda \left\{{\left(1+\epsilon z_3\right)}^m{e}^{\frac{z_3}{1+\epsilon z_3}}+\delta z_2^2\left(e^{-\alpha z_3}+2\gamma e^{-\alpha z_3}\right)\right\}\right)\end{array}\right)\\ \end{array}$$

(25)

with the set of initial conditions

$$\left(\begin{array}{l}{z}_1\left(0\right)\\ z_2\left(0\right)\\ z_3\left(0\right)\\ z_4\left(0\right)\end{array}\right)=\left(\begin{array}{l}0\\ D_1\\ \frac{z_4\left(0\right)e^{-\eta z_3}}{B{i}_1}\\ D_2\end{array}\right)$$

(26)

Since the derivatives are not specified at the initial point, $D_{1,2}$ are therefore given initial guess values that ensure the boundary conditions are okay at the wall y = 1. It is easy to see that the sufficient condition that confirms the existence of a unique solution to the system of Equations (24) and (25)is satisfied as long as the Lipchitz constant $\frac{\partial f_i}{\partial x_j}\le K$ is bounded. The first-order equations are then carefully coded and solved by Runge–Kutta using the Mathematica package.

5. Results and Discussion

In this section, a physical interpretation of the numerical results is presented and discussed as regards the flow and thermal behavior.

Table 1. Validation of numerical results for velocity $\delta =0.5,\eta =0.1,\alpha =0.3$, $\alpha =0.1,N=30,$ $\gamma =0.3,Gr=1,$ $B{i}_1=B{i}_2=10,$ $\lambda =0.4,\epsilon =0.1,m=0.5$.

$y$	$u{\left(y\right)}_{\mathbf{SCCM}}$	$u{\left(y\right)}_{\mathbf{RK}4}$	$\left|u{\left(y\right)}_{\mathbf{CWRM}}-u{\left(y\right)}_{\mathbf{RK}4}\right|$
0	−6.96195549259×10−18	0.000000000000000	6.961955492590054 × 10−18
0.1	0.04313941885340904	0.04313941993774605	1.084337016010739 × 10−9
0.2	0.07813634077297339	0.07813635190993772	1.11369643229775 × 10−8
0.3	0.10406977160813229	0.10406978921035206	1.760221977897824 × 10−9
0.4	0.12006990360382683	0.12006992567673559	2.207290876465872 × 10−8
0.5	0.12548527366302917	0.12548528888798596	1.522495679529001 × 10−9
0.6	0.12006990360382683	0.12006990533241987	1.728593046479432 × 10−9
0.7	0.1040697716081323	0.10406977022133693	1.386795372981808 × 10−9
0.8	0.07813634077297338	0.07813634228509685	1.512123468105919 × 10−9
0.9	0.043139418853409044	0.04313941955035403	6.96944987832459 × 10−10
1.0	−7.69196285588 × 10−18	−9.12736796 × 10−10	9.12736788740352 × 10−10

and

Table 2. Validation of numerical results for temperature.

$y$	$\theta {\left(y\right)}_{\mathbf{SCCM}}$	$\theta {\left(y\right)}_{\mathbf{RK}4}$	$\left|\theta {\left(y\right)}_{\mathbf{SCCM}}-\theta {\left(y\right)}_{\mathbf{RK}4}\right|$
0	0.022151849024484603	0.022151849374097372	3.496127692903528 × 10−10
0.1	0.04208489447368267	0.04208489042284373	4.050838935121259 × 10−9
0.2	0.05751544434774708	0.057515440006678895	4.341068185476082 × 10−9
0.3	0.06849268440471615	0.06849268031195316	4.092762989627019 × 10−9
0.4	0.07505930588225566	0.07505930197957113	3.902684536649659 × 10−9
0.5	0.07724467496042003	0.07724467142322092	3.537199116943057 × 10−9
0.6	0.07505930588225566	0.07505930282648429	3.055771372051374 × 10−9
0.7	0.06849268440471615	0.06849268167736071	2.727355438714163 × 10−9
0.8	0.05751544434774707	0.05751544180643439	2.541312688064678 × 10−9
0.9	0.04208489447368267	0.04208489220244641	2.271236268502896 × 10−9
1.0	0.022151849024484603	0.02215184702121924	2.003265362621187 × 10−9

represent the validation of the point collocation weighted residual method with the shooting-Runge–Kutta procedure for the dimensionless momentum and energy equations. A good agreement is observed between the two computations. This confirmed the uniqueness of the numerical solutions.

Table 3. Fast convergence of critical values by weighted residual method when $Gr=3,\delta =0.2,B{i}_1=10=B{i}_2,m=0.5,\epsilon =0.1$.

$N$	$\alpha$	$\sigma$	$\gamma$	$\eta$	${\lambda }_c$
5	0.1	0.1	0.3	0.1	3.369629994329486
10	0.1	0.1	0.3	0.1	3.329979429762467
15	0.1	0.1	0.3	0.1	3.3377915382719980
20	0.1	0.1	0.3	0.1	3.3375451276368135
25	0.1	0.1	0.3	0.1	3.3376391452777585
30	0.1	0.1	0.3	0.1	3.3376334297519272
35	0.1	0.1	0.3	0.1	3.3376356319779727

shows the convergence of the of the spectral collocation solution. The convergence of the solution to the highly nonlinear boundary value problem is seen to improve with an increasing number of terms $N$ used in the approximation.

Table 4. Thermal stability with values when $\delta =0.2,B{i}_1=10=B{i}_2,m=0.5,\epsilon =0.1$.

$\sigma$	$\alpha$	$\gamma$	$\eta$	$Gr$	${\lambda }_c$
0.1	0.1	0.3	0.1	3	3.337791538271998
0.3	0.1	0.3	0.1	3	3.294547099649542
0.5	0.1	0.3	0.1	3	3.251342585235665
0.1	0.3	0.3	0.1	3	3.327113416799457
0.1	0.5	0.3	0.1	3	3.3180031889548034
0.1	0.1	0.5	0.1	3	3.354693697420452
0.1	0.1	0.8	0.1	3	3.370237078964153
0.1	0.1	0.3	0.01	3	3.5150936154732557
0.1	0.1	0.3	0.001	3	3.5341579072708558
0.1	0.1	0.3	0.1	5	3.1911369048073452
0.1	0.1	0.3	0.1	7	3.0387490555917243

represents the computational result of the thermal stability analysis of the fluid. As observed from the plot, an increase in the nonlinear component of the Grashof number is seen to destabilize the flow thermal structure due increasing buoyancy force over viscous forces in the flow field. Additionally, an increase in the viscosity variation parameter is also observed to have a destabilizing effect on the fluid flow. This is true, since viscosity is expected to drop and fluid experiences shear thinning, which aids heat generation by kinetic theory of matter. Moreover, a rise in the third-grade parameter is seen to enhance thermal stability due to the thickening effect of the fluid. The thermal conductivity variation parameter increase is seen to discourage the thermal stability of the fluid flow. A similar behavior is seen with increasing values of the Grashof number due to increasing buoyancy forces.

Table 5. Validation with previous result when $B{i}_1=B{i}_2=10,Gr=0=\gamma =\eta ,\alpha =0.1$.

$y$	$u_1\left(y\right)$ Makinde & Aziz [34]	$u_2\left(y\right)$ Present Result	$\left|u_1\left(y\right)-u_2\left(y\right)\right|$
0	0.0000000000000000	1.9996146847 × 10−18	1.999614684734152 × 10−18
0.1	0.04500260897270728	0.04500261769036607	8.717658789292315 × 10−9
0.2	0.08000542199231503	0.08000543151384462	9.521529592548816 × 10−9
0.3	0.10500767322725557	0.10500768307611696	9.848861393102482 × 10−9
0.4	0.12000907000446212	0.12000908068447211	1.068000998749596 × 10−8
0.5	0.12500954265746725	0.12500954933196154	6.674494290592747 × 10−9
0.6	0.12000907650344211	0.12000908068447211	4.181029994443364 × 10−9
0.7	0.10500768080412677	0.10500768307611698	2.271990207081131 × 10−9
0.8	0.08000543108796186	0.08000543151384461	4.258827457359615 × 10−10
0.9	0.045002617803451495	0.04500261769036607	1.130854229702826 × 10−10
1.0	8.38288175864 × 10−10	2.13968421631 × 10−18	8.38288173724192 × 10−10

and

Table 6. Validation with previous result when $B{i}_1=B{i}_2=10,Gr=0=\gamma =\eta ,\alpha =0.1$.

$y$	${\theta }_1\left(y\right)-SRK4$ Makinde & Aziz [34]	${\theta }_2\left(y\right)-SCCM$ Present Result	$\left|{\theta }_1\left(y\right)-{\theta }_2\left(y\right)\right|$
0	0.0004166959732611298	0.00041669643079810236	4.57536972555978 × 10−10
0.1	0.0007242078568868499	0.0007242204230993113	1.256621246132922 × 10−8
0.2	0.0008700589609001264	0.0008700665516844415	7.59078431502528 × 10−9
0.3	0.0009242312498706879	0.0009242382030372055	6.953166517597101 × 10−9
0.4	0.0009367330606715411	0.0009367393712081855	6.310536644453524 × 10−9
0.5	0.0009375671371315745	0.0009375727826577018	5.645526127313893 × 10−9
0.6	0.000936734339945098	0.0009367393712081855	5.031263087466702 × 10−9
0.7	0.0009242337880723514	0.0009242382030372052	4.414964853774461 × 10−9
0.8	0.0008700627611325915	0.0008700665516844411	3.790551849616394 × 10−9
0.9	0.0007242172747296661	0.0007242204230993107	3.148369644634038 × 10−9
1.0	0.0004166939431660092	0.0004166964307981019	2.487632092722583 × 10−9

revealed the accuracy of the solution when compared with already-published results in the special case when the fluid is Newtonian, non-reactive, with constant physical properties.

Figure 2 represents the variation of third-grade material effect with velocity, temperature, entropy generation, and Bejan number, as shown in Figure 2a–d, respectively. The thickening effect is manifested in Figure 2a as the third-grade material parameter grows. Likewise in Figure 2b, thickening of the fluid is seen to make the non-Newtonian fluid colder. This is physically correct because the kinetic energy of the flow decreases as a result, the intermolecular bond is expected to increase. Therefore, the activation energy requirement will be much higher to induce chemical reaction. The result from Figure 2c reveals that as the third-grade material parameter rises, there is drop in the entropy generation at the cooled walls. This is reasonable, because a rise in the third-grade material parameter has been shown to reduce the flow velocity and temperature maximum in Figure 2a,b. Therefore, available energy for work is maximum with the thickened fluid. Finally, in Figure 2d the result of the heat irreversibility ratio credited Bejan is presented. If we suppose that the fluid velocity and density are both constant, then we have a case of slab with Be = 1. However, for a variable viscous third-grade flow, we have seen that velocity is a decreasing function of the third-grade parameter. As a result, the kinetic energy of the flow drops with increasing value of the parameter. Therefore, heat transfer irreversibility is expected to gain dominance over frictional irreversibility as shown in the plot.

Figure 3a–d represent the impact of buoyancy forces on the flow field on the velocity, temperature, entropy generation rate, and the Bejan number respectively. From Figure 3a, the velocity profile shows that when Gr = 1 implies a simple case when viscous force and buoyancy force are equal. Further increase in the Grashof number implies that the contribution of viscous force is weaker when compared with buoyancy, therefore, the flow velocity increases significantly in response to the nonlinear buoyancy forces. A similar behavior is seen in the profile for fluid temperature (see Figure 3b) resulting from increase in the kinetic energy of the fluid. Consequently in Figure 3c, the entropy profile is expected to increase due to increase in both flow and temperature maximum. In Figure 3d, rise in Grashof number implies increasing temperature as seen in Figure 3b. Therefore, increasing temperature translates to a decreased rate of heat transfer, this shows that irreversibility due to frictional interactions is more when compared with heat transfer irreversibility as revealed in Figure 3d. The graphical representations of the effect of the nonlinear convection component of the Grashof number on the flow velocity, temperature, entropy generation, and Bejan number respectively depicted in Figure 4a–d. The point when $\sigma =0$ corresponds to the linearized buoyancy parameter when the nonlinear effect is neglected. As shown in the Figure 4a, this parameter has a significant effect on the fluid flow and cannot be neglected. Also, in Figure 4b, an increased value of the parameter is also seen to elevate the fluid temperature maximum. Similar effect is noticed in Figure 4c, the rate of entropy generation improves at both convective walls while in Figure 4d, the frictional irreversibility decreases within the flow channel, due to a decrease in heat transfer.

The influence of the reaction kinetics coefficient ($\lambda$) on the flow and thermal structure is presented in Figure 5a. It is evidently shown that the velocity distribution enhanced for larger $\lambda$. This is due to the thinning effect on the fluid viscosity and the heat transfer from the strongly exothermic reaction to the fluid. In Figure 5b, similar behavior is experienced with the fluid temperature, except at the cooled walls with increasing flow and temperature. As expected the entropy generation rate is to rise for increased values of $\lambda$ as shown in Figure 5c, while larger $\lambda$ significantly increase the Bejan profile. This in turn encourages heat transfer irreversibility over fluid friction irreversibility as seen in Figure 5d.

Finally, Figure 6a–c shows the bifurcation plots for different reaction kinetics (m = −2, 0, 0.5). As observed from these plots, the blow-up point increases with increasing values of the exponent, thus thermal stability is enhanced.

6. Concluding Remarks

In this study, a time-independent flow and heat transfer of reactive third-order fluid has been investigated numerically, taking dependence of some physical properties on temperature into consideration. The nonlinear equations were solved numerically by using a spectral Chebyshev weighted residual method and validated by shooting-Runge–Kutta approach. The main contributions to knowledge from the study are listed as follows:

-

Flow, temperature, and entropy generation are enhanced with increasing values of the Grashof number, the quadratic component of buoyancy and Frank-Kameneskii parameter but reduces with increasing third-grade material parameter.

-

Increasing values of third-grade parameter encourages the thermal stability of the flow, while increasing values of the linear and nonlinear buoyancy parameter destabilizes the flow.

-

The Grashof number increase encourages the early occurrence of thermal runaway and exergy loss in the flow domain.

Future directions on this study are not limited to the developing flow, mass transfer, and elasto-hydrodynamics aspects of the work.