A Theoretical Analysis for Mixed Convection Flow of Maxwell Fluid between Two Infinite Isothermal Stretching Disks with Heat Source/Sink

Abstract

The aim of this current contribution is to examine the rheological significance of Maxwell fluid configured between two isothermal stretching disks. The energy equation is also extended by evaluating the heat source and sink features. The governing partial differential equations (PDEs) are converted into the ordinary differential equations (ODEs) by using appropriate variables. An analytically-based technique is adopted to compute the series solution of the dimensionless flow problem. The convergence of this series solution is carefully ensured. The physical interpretation of important physical parameters like the Hartmann number, Prandtl number, Archimedes number, Eckert number, heat source/sink parameter and the activation energy parameter are presented for velocity, pressure and temperature profiles. The numerical values of different involved parameters for skin friction coefficient and local Nusselt number are expressed in tabular and graphical forms. Moreover, the significance of an important parameter, namely Frank-Kamenetskii, is presented both in tabular and graphical form. This particular study reveals that both axial and radial velocity components decrease by increasing the Frank–Kamenetskii number and stretching the ratio parameter. The pressure distribution is enhanced with an increasing Frank–Kamenetskii number and stretching ratio parameter. It is also observed that thetemperature distribution increases with the increasing Hartmann number, Eckert number and Archimedes number.

1. Introduction

The mixed convection flow is the combination of both coupled free and forced convection, and is a topic of particular interest from an engineering (aerospace and chemical engineering) point of view in the past few years. A diverse significance of such a phenomenon may appear in various electronic devices, nuclear reactors, food industries, energy storage, era of astrophysics, lubrication phenomenon, fire control, chemical metallurgical, etc. The phenomenon of free convection is resulted due to the temperature difference in fluid particles associated with isothermal stretching disks.

The involvement of magnetic force in the heat transfer processes between stretching disks is termed as forced convection. In a mixed convection flow, the Archimedes number represents the comparative contribution of natural to forced convection. It is a well justified fact that the phenomenon of free convection becomes more prevailing over forced convection when the Archimedes number is larger than unity. In the modern era of science, the flows caused by heat supplied in the presence of transport processes which occurred due to chemical reactions gained the attention of investigators due to numerous applications in several industrial processes. Arrhenius kinetics is adopted for modeling such reactions, where the flow is thermally obsessed by exothermic surface reaction. Maleque [1] studied the effects of exothermic/endothermic chemical reactions in the presence of energy activation over a porous flat plate. The impact of nonlinear thermal radiation and activation energy in the flow of Cross nanofluid has been reported by Khan et al. [2]. Shafique et al. [3] examined the flow of Maxwell fluid along with activation energy features in a rotating frame. A numerically-based continuation for viscous fluid flow in the presence of activation energy and slip factors has been pointed out by Awad et al. [4]. Another interesting contribution on the flow of viscoelastic fluid in presence of activation energy was investigated by Hsiao [5]. According to this study, the obtained observations can be used to enhance the manufacturing and thermal extrusion systems. The mixed convection flow on chemically reactive surfaces for external flow in the presence of porous medium was investigated by Merkin and Mahmood [6]. Similar studies were also performed by Minto et al. [7] for a vertical surface. We also acknowledge the interesting study presented by Chou and Tsern [8], in which they presented experimentally-based results regarding mixed convection flow in a horizontal channel with constant heat flux conditions.

In recent years, the stretched flows of electrically conducting materials under the influence of magnetic force have attained attention due to diverse engineering and medical applications. Some valuable applications of this phenomenon may include nuclear reactors, fission and fusion reactions, plasma, metallurgical processes, the exploration of oil, thermal conductors, magnetohydrodynamic (MHD) generators, etc. The MHD flow passing in arteries is important because of diverse physiological processes. For example, the flow of blood can be effectively controlled via an addition of the mixing of samples, heat transportation and interaction of the magnetic field. Many authors performed an extensive analysis regarding the MHD flow of various fluid models with different geometries. Nadeem et al. [9] investigated the impact of magnetic force in viscous nanofluid flow configured by a curved surface. Ahmed et al. [10] performed some numerical computations while explaining the thermophysical consequences in nanofluid flow subjected to magnetic force. Khan et al. [11] successfully obtained the dual solution for the combined heat and mass flow of magnetized nanoparticles over a curved surface. The oscillatory flow of micropolar nanofluid subjected to magnetic force has been numerically inspected by Sadiq et al. [12].

The study of non-Newtonian fluids is important due to their wide range of applications in engineering, physiology, the chemical and petroleum industries. The non-Newtonian fluid models capture a nonlinear relationship between shear stress and deformation rate in contrast to the viscous materials. The traditional examples of such fluids include paints, blood, paste, jell, apple source, etc. The non-Newtonian boundary layer flow due to stretching surfaces has been paid a great attention by scientists due to interesting industrial and engineering applications like glass fiber manufacturing, paper production, plastic films, crystal growing and in the processing of cooling bath of metallic sheets. In order to study the physical properties of non-Newtonian fluids, various models have been introduced in the literature. The classification of non-Newtonian models can be referred as rate type, differential type and integral type fluids. In the category of rate type, Maxwell fluid is considered as a subclass of rate type liquids which accomplishes the relaxation time features. The examples of Maxwell fluid include crude oil, toluene, polymer solution, etc. Haris [13] suggested the boundary layer equations for two-dimensional flow of Maxwell fluid. After that the analysis of boundary layer flow and heat transfer over a stretching surface by using the constitutive equation of Maxwell fluid was carried out by several researchers. For instance, Hayat et al. [14] discussed the series solution of upper-convected Maxwell fluid over a porous stretching plate.

The effects of thermal radiation on the MHD flow of a Maxwell fluid over a stretching surface were examined by Aliakbar et al. [15]. Two-dimensional stagnation-point flow of upper-convected Maxwell fluid (UCM) over a stretching sheet has been determined by Hayat et al. [16]. They used the homotopy analysis method (HAM) to solve the resulting nonlinear differential equations. Prasad et al. [17] discussed the effects of temperature-dependent viscosity, thermal conductivity and internal heat generation/absorption features in the MHD flow of upper-convected Maxwell fluid configured by a stretched surface. Khan et al. [18] examined the flow of Maxwell fluid in a channel with oscillating walls under the action of a magnetic field. The analysis for Maxwell fluid in the presence of a heat transfer phenomenon over coaxially rotating disks has been depicted by Ahmed et al. [19].

The fluid flow encountered by a rotating and stretching disk has gained serious importance in the last years due to a large number of physical applications for both physical and theoretical aspects. Some emerging applications of such flows includes a rotor-stator system, MHD generators, turbine engines, aircraft engines, spin coating, centrifugal pumps, flow-through swept wings, shrouded-disks rotation, rotating electrodes, centrifuges, hydraulic press, boilers, condensers, etc. Merkin and Chaudhary [20] investigated the flow of viscous fluid induced by stretching a disk in the presence of an exothermic surface reaction. Gorder et al. [21] reported the analytical solution for flow encountered by stretching disks. Khan et al. [22] studied the mixed convection flow induced by exothermal and isothermal stretching disks analytically.

In the present analysis, we study an incompressible mixed convection flow of Maxwell fluid between infinite isothermal stretching disks in the presence of heat absorption/generation, activation energy and chemical reaction features. In fact, this work is the extension of Gorder et al. [21] in three directions: Firstly, by considering Maxwell fluid, secondly by including activation energy consequences, and lastly by taking heat source and sink features. Considering the literature survey, it is noted that this present analysis has not been investigated yet and presented for the first time in literature. The study of the mixed convection flow of non-Newtonian fluid encountered enormous applications in nuclear engineering, chemical engineering and petroleum industries. The considered flow problem contained the impact of activation energy, which includes diverse industrial and engineering significance, like oil emulsion, food processing, chemical processes and geothermal reservoirs. The problem is solved analytically via the homotopy analysis method, and the results are discussed through pictorial and tabular representations.

2. Mathematical Formulation

In the current analysis, a non-Newtonian fluid is configured between two infinite stretching disks. It is assumed that flow is axisymmetric and steady. The rheological aspects of non-Newtonian material have been deliberated by using the famous Maxwell fluid model which occupies the space $0<z<d$. The disks are separated distance d from each other as shown in Figure 1. The flow is generated due to the stretching of both disks in the radial direction. It is assumed that both (upper and lower) disks are isothermal in nature at temperatures $T_1$ and $T_2$, respectively. The analysis is performed by opting for cylindrical coordinates $\left(r,\theta ,z\right)$. All the involved expressions are independent of $\theta$ due to axisymmetry. Following Merkin et al. [20], the expressions for first order non-isothermal reaction are represented in following form

$$A\to B+heat, rate=k_0{a}_0{e}^{-E/R_{1}T}.$$

(1)

These above relations are known as Arrhenius kinetics, where E signifies the activation energy, B is a product species, $R_1$ is the gas constant, $k_0$ is the chemical reaction, $a_0$ the reactant concentration and T is the fluid temperature. The flow equations for the axisymmetric flow of Maxwell fluid can be expressed as [20,21,22]:

$$\frac{1}{r}\frac{\partial }{\partial r}\left(ru\right)+\frac{\partial w}{\partial z}=0,$$

(2)

$$\begin{array}{cc}\hfill u\frac{\partial u}{\partial r}+w\frac{\partial u}{\partial z}& =-\frac{1}{\rho }\frac{\partial p}{\partial r}+\nu \left(2\frac{{\partial }^{2}u}{\partial r^2}+\frac{{\partial }^{2}w}{\partial r\partial z}+\frac{{\partial }^{2}u}{\partial z^2}+\frac{2}{r}\frac{\partial u}{\partial r}-2\frac{u}{r^2}\right)\hfill \\ & -{\lambda }_1\left(w^2\frac{{\partial }^{2}u}{\partial z^2}+2uw\frac{{\partial }^{2}u}{\partial r\partial z}+u^2\frac{{\partial }^{2}u}{\partial r^2}\right)+\frac{\delta B_0{}^2}{\rho }\left(-u-{\lambda }_{1}w\frac{\partial u}{\partial z}\right)\hfill \\ & +g\beta \left[\left(T-T_0\right)+{\lambda }_1\left(u\frac{\partial T}{\partial r}+w\frac{\partial T}{\partial z}-\frac{\partial u}{\partial r}\left(T-T_0\right)\right)\right]\hfill \end{array}$$

(3)

$$\begin{gathered} u\frac{\partial w}{\partial r}+w\frac{\partial w}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial z}+\nu \left(\frac{{\partial }^{2}w}{\partial r^2}+\frac{{\partial }^{2}u}{\partial r\partial z}+2\frac{{\partial }^{2}w}{\partial z^2}+\frac{1}{r}\frac{\partial w}{\partial r}+\frac{1}{r}\frac{\partial u}{\partial r}\right) \\ -{\lambda }_1\left(w^2\frac{{\partial }^{2}u}{\partial z^2}+2uw\frac{{\partial }^{2}u}{\partial r\partial z}+u^2\frac{{\partial }^{2}u}{\partial r^2}\right), \end{gathered}$$

(4)

$$\begin{array}{cc}\hfill \rho c_p\left(u\frac{\partial T}{\partial r}+w\frac{\partial T}{\partial z}\right)& \hfill \\ & =2\mu ({\left(\frac{\partial u}{\partial r}\right)}^2+{\left(\frac{\partial u}{\partial z}\right)}^2+2\left(\frac{\partial w}{\partial r}\right)\left(\frac{\partial u}{\partial z}\right)+{\left(\frac{u}{r}\right)}^2+{\left(\frac{\partial w}{\partial r}\right)}^2+{\left(\frac{\partial w}{\partial z}\right)}^2)\hfill \\ & + K_T\left(\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{{\partial }^{2}T}{\partial z^2}\right)-\frac{\delta B_0{}^{2}u}{\rho }+Q_0\left(T-T_0\right)+Q{k}_0{a}_0{e}^{-E/R_{1}T},\hfill \end{array}$$

(5)

in which $u$ and $w$ are velocity components in the $r$ and $z$ directions, ${\lambda }_1$ is the relaxation time,$p$ is the fluid pressure, $\rho$ is the characteristic density, $\nu$ is the kinematic viscosity, $K_T$ is the thermal conductivity of the fluid, $T_0$ the reference temperature given by $T_0=\frac{T_1+T_2}{2}$, β denotes the thermal expansion coefficient, $Q_0$ denotes the heat generation/absorption coefficient while $Q$ stands for the exothermicity factor. The imposed boundary conditions associated with the current flow problem are:

$$\begin{gathered} u=ar, w=0, p=\frac{a\mu \beta r^2}{4d^2}, at z=0 \\ u=cr, w=0, p=0, at z=d, \\ T=T_1 at z=0, T=T_2, at z=d. \end{gathered}$$

(6)

In order to obtain the dimensionless form of above equations, we introduce the following similarity variables [21,22]:

$$\begin{gathered} u=arF\left(\eta \right), w=adH\left(\eta \right), \eta =\frac{z}{d} \\ u=arF\left(\eta \right), w=adH\left(\eta \right), \eta =\frac{z}{d} \\ p=a\mu \left(P\left(\eta \right)+\frac{{\beta }_1{r}^2}{4d^2}\right), \theta \left(\eta \right)=\frac{E\left(T-T_0\right)}{R_1{T}_0{}^2}. \end{gathered}$$

(7)

The above transformations lead to the following system:

$$\begin{array}{cc}\hfill H^{‴}+\beta +\frac{\mathrm{R}}{2}({H^{\prime }}^2& -2H{H}^{″})\hfill \\ & ={\lambda }_{1}aR\left(H^2{H}^{‴}-H{H}^{\prime }{H}^{″}\right)-M\left(H^{\prime }+{\lambda }_{1}H{H}^{″}\right)+A_r[2\theta \hfill \\ & + {\lambda }_{1}a\left(2H{\theta }^{\prime }+H^{\prime }\theta \right)],\hfill \end{array}$$

(8)

$${\theta }^{″}-RPrH{\theta }^{\prime }+\frac{EcPr}{ϵ}\left(\frac{{\delta }^2}{2}{\left(H^{″}\right)}^2+3{H^{\prime }}^2\right)+\alpha \theta +K\exp \left(\frac{-1}{ϵ\left(1+ϵ\theta \right)}\right)=0,$$

(9)

$$P^{\prime }=H^{″}-RH{H}^{\prime }-\lambda R{H}^2{H}^{″},$$

(10)

$$\begin{array}{c}H\left(0\right)=0, H\left(1\right)=0, H^{\prime }\left(0\right)=-2, H^{\prime }\left(1\right)=-2\gamma ,\\ \theta \left(0\right)=R_T , \theta \left(1\right)=-R_T, P\left(0\right)=0. \end{array}\}$$

(11)

where the stretching rate constant is $\gamma$, the Reynolds number $R$, the Hartmann number is $M$, Grashoff number $G_r$, heat source/sink parameter $\alpha ,$ the Prandtl number $Pr$, Eckert number $Ec$, the Frank–Kamenetskii number $K,$ constant temperature parameter $R_T$, activation energy parameter $ϵ$, the Archimedes number $A_r$ and the dimensionless distance $\delta$ are defined as:

$$\begin{array}{c}\gamma =\frac{c}{a},R=\frac{a{d}^2}{v}, M= \frac{\sigma B_0{}^2}{a\rho }, G_r= \frac{g\beta T_0{r}^3}{{\nu }^2}, \alpha =\frac{Q_0{d}^2}{K_T} , Pr=\frac{\mu c_p}{K_T}, \\ Ec=\frac{a{d}^2}{C_p{T}_0},K=\frac{Q{k}_0{a}_0{d}^2}{K_TϵT_0}{e}^{-\frac{1}{ϵ}}, R_T=\frac{T_1-T_0}{ϵT_0}, ϵ= \frac{R_1{T}_0}{E},A_r= \frac{G_r}{R^2},\delta =\frac{r}{d}\end{array}\}$$

(12)

By differentiating Equation (8) with respect to similarity variable $\eta$, we have

$$\begin{array}{cc}\hfill H^{\left(iv\right)}-RH{H}^{‴}& =\lambda R\left(H{H}^{\prime }{H}^{‴}+H^2{H}^{\left(iv\right)}-{H^{\prime }}^2{H}^{″}-H{H^{″}}^2\right)\hfill \\ & -MR \left[H^{″}+\lambda \left(H^{\prime }{H}^{″}+H{H}^{‴}\right)\right]\hfill \\ & +\frac{R{A}_rϵ}{{\delta }^4}\left[2{\theta }^{\prime }+\lambda \left(3H^{\prime }{\theta }^{\prime }+2H{\theta }^{″}+H^{″}\theta \right)\right]. \hfill \end{array}$$

(13)

where $\lambda ={\lambda }_{1}a$, is the Deborah number. First of all, we solve Equation (12) subject to the boundary conditions (11) and then $\beta$ can be evaluated by using Equation (8).

2.1. Skin Friction Coefficient

The expression for shear stress ${\tau }_w$ on the surface of the stretching disk is defined as [21,22]:

$${\tau }_w={\tau }_{rz}|{ }_{z=0},$$

(14)

The skin friction coefficients $R{C}_{1f}$ and $R{C}_{2f}$ at the lower and upper disks are:

$$C_{1f}=\frac{{\tau }_w}{\frac{1}{2}\rho {\left(\delta r\right)}^2}=\frac{{\tau }_{rz}|{ }_{z=0}}{\frac{1}{2}\rho {\left(\delta r\right)}^2}=-R^{-1}{H}^{″}\left(0\right),$$

(15)

$$C_{2f}=\frac{{\tau }_w}{\frac{1}{2}\rho {\left(\delta r\right)}^2}=\frac{{\tau }_{rz}|{ }_{z=d}}{\frac{1}{2}\rho {\left(\delta r\right)}^2}=-R^{-1}{H}^{″}\left(1\right).$$

(16)

2.2. Local Nusselt Number

The mathematical expressions for the local Nusselt number is represented as:

$$q_w=-\left(K_T\frac{\partial T}{\partial z}\right)=-\frac{K_T{R}_1{T}_0^2}{Ed}{\theta }^{\prime }\left(\eta \right),$$

(17)

The dimensionless form of this local Nusselt number at both (lower and upper) disks is [21,22]:

$$N_{1u}=\frac{d{q}_w}{K_TϵT_0}=-\frac{d{K}_T\frac{\partial T}{\partial z}|{ }_{z=0}}{K_TϵT_0}=-\frac{d{K}_T{R}_1{T}_0^2}{Ed{K}_TϵT_0}{\theta }^{\prime }\left(0\right)=-{\theta }^{\prime }\left(0\right),$$

(18)

$$N_{2u}=\frac{d{q}_w}{K_TϵT_0}=-\frac{d{K}_T\frac{\partial T}{\partial z}|{ }_{z=d}}{K_TϵT_0}=-\frac{d{K}_T{R}_1{T}_0^2}{Ed{K}_TϵT_0}{\theta }^{\prime }\left(1\right)=-{\theta }^{\prime }\left(1\right).$$

(19)

3. HomotopyAnalysis Method

In our modern era of scientific research, many physical and engineering problems are modeled in the form of highly nonlinear differential equations which always remain challenging for mathematicians to suggest the analytical or numerical solutions. Among different analytical techniques, thehomotopy analysis method is one which can be used to compute the analytic solution of such problems with excellent convergence. This technique is free of a complicated discretization procedure like numerical methods. This analytical technique is free of any small or large parameter constraints. This method was originally introduced by Liao [23], and later on many researchers used this method for their solutions of various problems [24,25,26,27,28,29,30,31,32,33,34,35,36]. The initial guesses for $H\left(\eta \right)$ and $\theta \left(\eta \right)$ are given by:

$$H_0\left(\eta \right)=2\eta \left(1-\eta \right)\left(\left(1+\gamma \right)\eta -1\right),$$

(20)

$${\theta }_0\left(\eta \right)=\eta ,$$

(21)

Defining auxiliary linear operators

$$L_H\left[y\right]=\frac{d^{4}y}{d{\eta }^4},$$

(22)

$$L_{\theta }\left[y\right]=\frac{d^{2}y}{d{\eta }^2},$$

(23)

Satisfying

$$L_H\left[C_1+C_2\eta +C_3{\eta }^2+C_4{\eta }^3\right]=0,$$

(24)

$$L_{\theta }\left[C_5+C_6\eta \right]=0,$$

(25)

where $C_i\left(i=1-6\right)$ are constants.

4. Convergence of Obtained Solution

It is a well-established fact that the convergence rate of HAM solutions is strictly based on non-zero auxiliary parameters ${\hslash }_H$ and ${\hslash }_{\theta }$. The suitable selection of these parameters is quite useful for adjusting and controlling the obtained solution. The admissible range of these auxiliary parameters, the $\hslash -$curves for velocity and temperature distributions, is displayed in Figure 2a,b. These figures clearly demonstrate that the suitable values of ${\hslash }_H$ and ${\hslash }_{\theta }$ can be selected from $-2.0\le {\hslash }_H\le -0.1$ and $-1.8\le {\hslash }_{\theta }\le -0.3.$ For present computations, the optimal values of ${\hslash }_H$ and ${\hslash }_{\theta }$ are taken ${\hslash }_H=-1$ and ${\hslash }_{\theta }=-1.08.$ The accuracy of obtained solution against various values of emerging parameters is shown in

Table 1. The convergence analysis of the homotopic solution with $R=5, \gamma =0.5, M=0.5, \lambda =0.2, A_r=2, \delta =3, Ec=1, Pr=1, \alpha =0.5, K=0.01, R_T=2, ϵ=0.5, {\hslash }_H=-1$ and ${\hslash }_{\theta }=-1.08.$

Order of Approximation	${\mathit{H}}^{″}\left(0\right)$	${\mathit{\theta }}^{\prime }\left(0\right)$
11	9.79594	−1.91738
14	9.79619	−1.91686
16	9.79643	−1.91634
18	9.79667	−1.91582
20	9.79717	−1.91461
25	9.79717	−1.91461
30	9.79717	−1.91461

. It is seen that the accuracy of the HAM solution is obtained at the twentieth order of approximations.

5. Validation of Solution

Before performing detailed graphical computations for flow parameters, we first compare our results with Gorder et al. [21] as a limiting case in

Table 2. Comparison of $H\left(\eta \right) \mathrm{and} H^{\prime }\left(\eta \right)$ for different values of $\eta$ with $\lambda =0, R=0$ and $M=0.$

$\mathit{\eta }$	Gorder et al. [21]		Present Result (HAM)
	$\mathit{H}\left(\mathit{\eta }\right)$	${\mathit{H}}^{\prime }\left(\mathit{\eta }\right)$	$\mathit{H}\left(\mathit{\eta }\right)$	${\mathit{H}}^{\prime }\left(\mathit{\eta }\right)$
0.0	0.000	−2.00	0.000	−2.00
0.2	−0.224	−0.360	−0.224	−0.360
0.4	−0.192	0.560	−0.192	0.560
0.6	−0.048	0.760	−0.048	0.760
0.8	0.064	0.240	0.064	0.240
1.0	0.000	−1.00	0.000	−1.000

. It is noted that an excellent accuracy of our results has been noted with these reported studies. Also Figure 3 shows the comparison of present results for the velocity profile computed via the homotopy analysis method for various values of the stretching ratio parameter with Gorder et al. [21]. It is found that present results have shown a convincible accuracy with Gorder et al. [21].

6. Results and Discussion

The formulated ordinary differential equations are targeted analytically via the homotopy analysis scheme. The aim of this section is to examine the physical significance of each physical parameter on velocity, pressure and temperature distributions.

6.1. Velocity Distribution

Figure 4, Figure 5 and Figure 6 are plotted to capture the influence of various parameters like the stretching ratio $\gamma$, Deborah number $\lambda$, Reynolds number $R$, Prandtl number $Pr$, heat source/sink parameter $\alpha$, constant temperature parameter $R_T,$ the Eckert number Ec, activation energy parameter $ϵ,$ Frank–Kamenetskii number $K,$ Hartmann number $M$ and the Archimedes number $A_r$ on the $r-$velocity component $H^{\prime }\left(\eta \right)$ and $z-$direction velocity component $H\left(\eta \right).$Figure 4a–f presents the effect of the Deborah number $\lambda ,$ stretching ratio $\gamma ,$ Eckert number $Ec,$ Prandtl number $Pr,$ dimensionless distance parameter $\delta$ and the Frank–Kamenetskii number $K$ on ther- and z-components of velocity. Figure 3a prescribed the outcomes of the Deborah number $\lambda$ on thez component of velocity. It is observed that the r-component of velocity increases up to a certain range and later on decreases slightly. It can be justified physically, as the Deborah number is directly proportional to the relaxation time. In fact, it is the associated with the fluid relaxation time to the observation time. The smaller values of the Deborah number represent the viscous nature of fluid while a material having a higher Deborah number represents the solid nature of fluid. The effects of the stretching ratio $\gamma$, Eckert number $Ec,$ dimensionless distance $\delta$ and the Frank–Kamenetskii number $K$ on $H\left(\eta \right)$ and $H^{\prime }\left(\eta \right)$ has been expressed in Figure 4b–e. From all these figures, it is noted that both $H\left(\eta \right)$ and $H^{\prime }\left(\eta \right)$ are enhanced by varying these parameters. Figure 4f manifested the influence of the Hartmann number on the z component of velocity. A decay in the z component of velocity is observed for intensifying values of the Hartmann number. Physically, the larger values of this Hartmann number attributed strong drag force which resists the amplitude of flow.Figure 4g–h determined the effects of the Archimedes number $A_r$ and Reynolds number $R$ on axial and radial velocities. A retarded distribution of both components has been resulted with the variation of all these parameters. Since the Reynolds number represents the ratio of inertial force to viscous, therefore higher values of $R$ become associated with larger inertial force which decay the velocity distribution effectively.

6.2. Pressure Distribution

The variation in pressure distribution $P\left(\eta \right)$ for various values of the stretching ratio parameter $\gamma ,$ dimensionless distance $\delta ,$ Reynolds number $R,$ Archimedes number $A_r,$ constant temperature parameter $R_T,$ activation energy parameter $ϵ,$ Deborah number $\lambda$, Hartmann number $M$ and the Frank–Kamenetskii number $K$ is discussed in Figure 5a–h. Figure 5a,b show the change in $P\left(\eta \right)$ for diverse values of $\gamma$ and $\delta .$ It is noted that pressure is an increasing function of $\gamma$ and $\delta$ up to a specific height, and later on decreases gradually. However, a decreasing trend has been observed for maximum values of the Reynolds number $R$ and the Archimedes number $A_r$ (Figure 5c,d). The graphical explanation for the Deborah number $\lambda$ and the Hartmann number $M$ is presented in Figure 5e,f. The Deborah number specified the relaxation time to the observation time ratio which means that maximum values of $\lambda$ correspond to larger relaxation time due to which the pressure distribution declined. Similarly, a decreasing trend in the pressure distribution is due to the fact that the Hartmann number is associated with Lorentz force, which efficiently controls the pressure distribution in the whole domain. Figure 5g determines the influence of the Frank–Kamenetskii number $K$ on pressure distribution $P\left(\eta \right).$ A retarded pressure distribution has been examined with the variation of $K.$ With the increase of $K,$ the pressure distribution decreases up to maximum level.

6.3. Temperature Distribution

In order to examine the impact of the Hartmann number $M,$ heat source/sink parameter $\alpha$, Eckert number $Ec,$ stretching ratio parameter $\gamma ,$ Archimedes number $A_r,$ dimensionless distance $\delta$ and the Reynolds number $R$ on temperature distribution $\theta \left(\eta \right),$Figure 6a–g are prepared. Figure 6a captured the consequences of the Hartmann number $M$ on temperature distribution $\theta \left(\eta \right).$ As expected, an enhanced temperature distribution is observed for larger values of $M$ due to the interaction of Lorentz force. From Figure 6b, again an increment in temperature distribution has been noted for maximum values of the Eckert number $Ec.$ The physical consequences of such trend may be attributed as heat due to viscous dissipation of fluid enhanced, due to which results an increment in $\theta \left(\eta \right)$. Figure 6c,d portrayed the impact Archimedes number $A_r$ and activation energy parameter $ϵ$ on $\theta \left(\eta \right).$ It is seen that temperature distribution enlarges with increasing both parameters. The activation energy plays a significant role in enhancement of many reaction processes. Figure 6e reports the influence of the heat source/sink constant on $\theta \left(\eta \right).$ It is noted that $\theta \left(\eta \right)$ increases in the case of heat source case $(\alpha >0)$, while the opposite trend is noted for the heat sink case $(\alpha <0)$. The physical aspect of such a trend may attribute, as in the case of heat source, more heat is added to the system, due to which the temperature distribution improved. On the contrary, due to the heat sink, heat is removed from the whole system which turns down the temperature distribution efficiently. From Figure 6f,g, a declining temperature distribution has been observed with maximum variation of dimensionless distance $\delta$ and Reynolds number $R.$

6.4. Physical Quantities of Interests

Figure 7a–c show the effect of different parameters like the stretching ratio parameter $\gamma ,$ Deborah number $\lambda$ and Hartmann number $M$ on the skin friction coefficient at upper and lower disks. From Figure 7a, a decreasing variation in this skin friction coefficient is examined with increasing $\gamma .$ On contrary, the skin friction coefficient at both level disks is increased for maximum values of the Deborah number $\lambda$(Figure 7b). Figure 7c reveals that the wall shear stress for different values of the Hartmann number $M$ is maximum at the upper level of the disk as compared to the lower level.

Figure 8a,b show the effects of Hartmann number $M$ and Eckert number $Ec$ on the local Nusselt number at lower and upper disks. The variation in local Nusselt number at the upper disk is larger for both parameters.

The numerical iteration in the wall shear stress at the upper level of the disk $R{C}_{1f}$ and lower level $R{C}_{2f}$ are discussed in

Table 3. Numerical variation in wall shear stress at both surfaces of moving disk.

$\mathit{\gamma }$	$\mathit{M}$	$\mathit{R}$	$\mathit{P}\mathit{r}$	$\mathit{E}\mathit{c}$	$\mathit{\lambda }$	$\mathit{\delta }$	$\mathit{\alpha }$	$\mathit{ϵ}$	${\mathit{A}}_{\mathit{r}}$	${\mathit{R}}_{\mathit{T}}$	$\mathit{K}$	Lower Disk	Upper Disk
0.2	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−116.204	−252.335
0.4	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−129.404	−262.989
0.6	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−142.249	−272.563
0.5	1.0	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−135.869	−267.915
0.5	1.5	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−149.413	−273.824
0.5	2.0	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−162.905	−281.867
0.5	01	1.0	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−1.63619	−38.7856
0.5	01	2.0	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−22.2454	−67.5845
0.5	01	3.0	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−44.9617	−108.970
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−135.869	−267.915
0.5	01	05	2.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−102.816	−437.258
0.5	01	05	3.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−69.7624	−606.601
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−135.869	−267.915
0.5	01	05	1.0	1.5	0.2	0.5	0.5	0.5	02	2.0	0.1	115.683	−430.368
0.5	01	05	1.0	2.0	0.2	0.5	0.5	0.5	02	2.0	0.1	367.236	−592.821
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−135.869	−267.915
0.5	01	05	1.0	1.0	0.4	0.5	0.5	0.5	02	2.0	0.1	−195.891	−369.982
0.5	01	05	1.0	1.0	0.6	0.5	0.5	0.5	02	2.0	0.1	−255.476	−472.018
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−135.869	−267.915
0.5	01	05	1.0	1.0	0.2	0.6	0.5	0.5	02	2.0	0.1	−62.2384	−77.8090
0.5	01	05	1.0	1.0	0.2	0.7	0.5	0.5	02	2.0	0.1	−23.2608	−44.7235
0.5	01	05	1.0	1.0	0.2	0.5	0.1	0.5	02	2.0	0.1	−130.344	−270.583
0.5	01	05	1.0	1.0	0.2	0.5	0.3	0.5	02	2.0	0.1	−133.108	−269.244
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−135.869	−267.915
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.1	02	2.0	0.1	−14.4518	−37.1367
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.3	02	2.0	0.1	−52.1849	−61.3395
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−89.1381	−99.7870
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	3.5	2.0	0.1	24.7854	−1205.14
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	4.0	2.0	0.1	195.575	−1786.25
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	4.5	2.0	0.1	445.195	−2541.17
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−135.869	−267.915
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.3	0.1	−185.374	−345.966
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.6	0.1	−247.552	−440.297
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−155.783	−257.838
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.2	−154.907	−257.367
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.3	−154.003	−256.903

. The wall shear stress gets minimum values for stretching rate constant $\gamma$, Reynolds number $R$, and Hartmann number $M$. It is noted that rate of wall shear stress is relatively slower at the lower portion of the disk for all parameters. The variation for various parameters on the local Nusselt number is portrayed in

Table 4. Numerical variation in local Nusselt number at both surfaces of moving disk.

$\mathit{\gamma }$	$\mathit{M}$	$\mathit{R}$	$\mathit{P}\mathit{r}$	$\mathit{E}\mathit{c}$	$\mathit{\lambda }$	$\mathit{\delta }$	$\mathit{\alpha }$	$\mathit{ϵ}$	${\mathit{A}}_{\mathit{r}}$	${\mathit{R}}_{\mathit{T}}$	$\mathit{K}$	Lower Disk	Upper Disk
0.2	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−43.2173	62.2902
0.4	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−43.1456	64.2379
0.6	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−43.1225	66.2978
0.5	1.0	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−43.1714	65.2526
0.5	1.5	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−49.2164	71.3174
0.5	2.0	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−55.2913	77.3855
0.5	01	2.0	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−1.93260	10.9702
0.5	01	3.0	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−9.76591	20.5937
0.5	01	4.0	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−20.8434	35.8615
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−34.5174	58.1590
0.5	01	05	2.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−59.2591	126.626
0.5	01	05	3.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−70.4176	209.251
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−34.5174	58.1590
0.5	01	05	1.0	1.5	0.2	0.5	0.5	0.5	02	2.0	0.1	−65.0929	88.3373
0.5	01	05	1.0	2.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−95.5317	118.241
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−37.1562	59.1911
0.5	01	05	1.0	1.0	0.4	0.5	0.5	0.5	02	2.0	0.1	−35.3967	58.5028
0.5	01	05	1.0	1.0	0.6	0.5	0.5	0.5	02	2.0	0.1	−33.6384	57.8153
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−34.5174	58.1590
0.5	01	05	1.0	1.0	0.2	0.6	0.5	0.5	02	2.0	0.1	−4.57816	17.6285
0.5	01	05	1.0	1.0	0.2	0.7	0.5	0.5	02	2.0	0.1	−2.16759	8.23666
0.5	01	05	1.0	1.0	0.2	0.5	0.1	0.5	02	2.0	0.1	−32.5620	56.3592
0.5	01	05	1.0	1.0	0.2	0.5	0.3	0.5	02	2.0	0.1	−33.5476	57.2673
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−34.5174	58.1590
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.1	02	1.0	0.1	−4.03154	4.78376
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.3	02	1.0	0.1	−7.31547	21.5100
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	1.0	0.1	−34.5174	58.1590
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	3.5	1.0	0.1	−128.255	180.574
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	4.0	1.0	0.1	−171.887	236.943
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	4.5	1.0	0.1	−221.713	301.095
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.0	0.1	−34.5174	58.1590
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.3	0.1	−43.9555	76.6222
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	2.6	0.1	−53.7965	97.8575
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	1.0	0.1	−35.1337	58.7998
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	1.0	0.2	−35.8183	59.5099
0.5	01	05	1.0	1.0	0.2	0.5	0.5	0.5	02	1.0	0.3	−36.5027	60.2180

. Again, the continuations are performed at both surfaces (upper surface $N_{1u}$ and lower surface $N_{2u}$). This physical quantity increases with $Pr$ and $Ec.$ Finally, the numerical values of Frank–Kamenetskii against different values of $\gamma , M, \lambda , R, Pr, Ec, \delta , \alpha , ϵ, A_r \mathrm{and} R_T$ is shown in

Table 5. The numerical values of Critical Frank–Kamenetskii number for various parameters.

$\mathit{\gamma }$	$\mathit{M}$	$\mathit{R}$	$\mathit{P}\mathit{r}$	$\mathit{E}\mathit{c}$	$\mathit{\delta }$	$\mathit{\alpha }$	$\mathit{ϵ}$	${\mathit{A}}_{\mathit{r}}$	${\mathit{R}}_{\mathit{T}}$	$\mathit{\lambda }$	${\mathit{K}}_{\mathit{c}}$
0.3	01	05	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.8144
0.5	01	05	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.8151
0.7	01	05	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.8170
0.5	01	05	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.8151
0.5	03	05	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.7897
0.5	05	05	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.7400
0.5	01	01	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.6770
0.5	01	02	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.6946
0.5	01	03	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.7235
0.5	01	05	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.8151
0.5	01	05	02	01	0.5	0.5	0.01	02	2.0	1.2	−10.9559
0.5	01	05	03	01	0.5	0.5	0.01	02	2.0	1.2	−11.0960
0.5	01	05	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.8151
0.5	01	05	01	02	0.5	0.5	0.01	02	2.0	1.2	−10.8064
0.5	01	05	01	03	0.5	0.5	0.01	02	2.0	1.2	−10.7978
0.5	01	05	01	01	2.0	0.5	0.01	02	2.0	1.2	−10.5643
0.5	01	05	01	01	2.5	0.5	0.01	02	2.0	1.2	−10.4735
0.5	01	05	01	01	3.0	0.5	0.01	02	2.0	1.2	−10.3631
0.5	01	05	01	01	0.5	0.3	0.01	02	2.0	1.2	−10.6729
0.5	01	05	01	01	0.5	0.5	0.01	02	2.0	1.2	−10.8151
0.5	01	05	01	01	0.5	0.7	0.01	02	2.0	1.2	−10.9575
0.5	01	05	01	01	0.5	0.5	0.1	02	2.0	1.2	−12.1179
0.5	01	05	01	01	0.5	0.5	0.2	02	2.0	1.2	−14.0721
0.5	01	05	01	01	0.5	0.5	0.3	02	2.0	1.2	−16.1551
0.5	01	05	01	01	0.5	0.5	0.01	10	2.0	1.2	−11.3132
0.5	01	05	01	01	0.5	0.5	0.01	50	2.0	1.2	−15.1132
0.5	01	05	01	01	0.5	0.5	0.01	100	2.0	1.2	−20.1335
0.5	01	05	01	01	0.5	0.5	0.01	02	01	1.2	−10.4104
0.5	01	05	01	01	0.5	0.5	0.01	02	02	1.2	−10.3151
0.5	01	05	01	01	0.5	0.5	0.01	02	03	1.2	−11.2584
0.5	01	05	01	01	0.5	0.5	0.01	02	2.0	0.5	−10.8071
0.5	01	05	01	01	0.5	0.5	0.01	02	2.0	1.0	−10.8127
0.5	01	05	01	01	0.5	0.5	0.01	02	2.0	1.5	−10.8189

. The variation in Frank–Kamenetskii constant is slower for $\lambda$ and $A_r.$

7. Conclusions

The axisymmetric flow of Maxwell fluid between two isothermal stretching disks is discussed in presence of source/sink and activation energy features. The mixed convection effects are implemented in the momentum equation. Analytical results are discussed by using the homotopy analysis method. The following observations are furnished:

• The wall shear stress decreases by increasing stretching parameter, Hartmann number, Reynolds number, Deborah number, activation energy parameter and constant temperature parameter. It means that tangential stresses increase by increasing stretching the ratio parameter, Hartmann number and Reynolds number. While the behavior of dimensionless distance and Frank–Kamenetskii number are quite the opposite.
• The pressure distribution is increased with variation of theFrank–Kamenetskii number and stretching ratio parameter.
• When the Deborah number $\lambda$ and Hartmann number increases, the wall shear stress at the lower disk increases while an opposite trend is found at the upper disk.
• It is observed that the surface heat transfer increases by increasing the stretching parameter and heat source/sink parameter.
• The rate of heat transfer decreases at the lower disk and increases at the upper disk by increasing the Hartmann number, Reynolds number, Archimedes number and activation energy parameter.