Entropy Analysis of 3D Non-Newtonian MHD Nanofluid Flow with Nonlinear Thermal Radiation Past over Exponential Stretched Surface

Abstract

The present study characterizes the flow of three-dimensional viscoelastic magnetohydrodynamic (MHD) nanofluids flow with entropy generation analysis past an exponentially permeable stretched surface with simultaneous impacts of chemical reaction and heat generation/absorption. The analysis was conducted with additional effects nonlinear thermal radiation and convective heat and mass boundary conditions. Apposite transformations were considered to transform the presented mathematical model to a system of differential equations. Analytical solutions of the proposed model were developed via a well-known homotopy analysis scheme. The numerically calculated values of the dimensionless drag coefficient, local Nusselt number, and mass transfer Nusselt number are presented, with physical insights. The graphs depicting the consequences of numerous parameters on involved distributions with requisite deliberations were also a part of this model. It is seen that the Bejan number is an increasing function of the thermal radiation parameter.

1. Introduction

Rapid growth in the usage of nanofluids in varied engineering arenas such as cancer therapy, finer coolants in nuclear reactors and computers, numerous electronic devices in military sectors [1], oil and water [2,3], rapid spry cooling, the food industry, vehicles and transformers, polymer extrusion, safe surgeries, quenching in foundries, and glass blowing [4] have encouraged scientists and researchers to scrutinize the numerous aspects of nanofluid flow past different geometries. Nanofluids are an amalgamation of nanoparticles of size < 100 nm (i.e., metals, nitrides, carbides, and nanotubes (single or multiwalled)) and orthodox fluids such as water, toluene, engine, kerosene oil, ethylene and triethylene glycol. Nanofluids are considered to be the best coolants in all engineering applications. It was Choi [5] who experimentally revealed that nanofluids possess enhanced thermal conductivity. Afterward, Buongiorno [6] framed a model pointing out that the heat transfer process is triggered in the case of nanofluids due to thermophoretic diffusion and the Brownian motion of nanoparticles. Later, Khan and Pop [7] discussed nanofluid flow past a stretched surface. Makinde and Aziz supported this model through convective boundary conditions over a stretched surface [8]. The impact of radiative nanoparticle flow with Lorentz force past a spongy semi annulus region was investigated by Sheikholeslami et al. [9]. Li et al. [10] examined the effect of adding nanoparticles in the process of solidification of Nano-enhanced phase change material (NEPCM) in the presence of thermal radiation. The flow of water-based nanofluid containing carbon nanotubes past a permeable medium with the effects of Darcy–Forchheimer was deliberated by Muhammad et al. [11]. Lu et al. [12] pondered on the flow of micropolar nanofluids in the presence of magnetohydrodynamic (MHD), homogeneous-heterogeneous reactions, mixed convection, and thermal radiation past a nonlinear stretched surface, numerically. The combined impacts of Arrhenius activation energy with heat and mass stratification on micropolar nanofluid flow in the presence of binary chemical reactions was studied by Ramzan et al. [13]. Lu et al. [14] did a tremendous study considering the three-layer vertical nanofluid model following a Buongiorno scheme. Li et al. [15], using Finite Element method (FEM), calculated the influence of the transportation of nanoparticles in a permeable enclosure due to electric force. Recently, some studies have discussed the applications of nanofluids in varied fields such as electrokinetic transport [16], heat transfer [17], and mass transfer [18].

The inevitable and dynamic applications of MHD nanofluids (e.g., wound treatment, gastric medications, targeted drug release, sterilized devices, magnetic resonance imagining (MRI), asthma treatment, and removal of tumors with hyperthermia) have attracted the attention of scientists and researchers in studying the related field of nanofluids. Sheikholeslami [19] deliberated the flow of copper oxide-water nanofluids in a spongy channel due to magnetic field effects using the Lattice Boltzmann method (LBM) method. In a recent article, Lu et al. [20] envisaged the model of MHD Carreau nanofluid flow over a radially stretching surface with allied impacts of nonlinear thermal radiation with newly established zero mass flux boundary conditions. They solved the proposed model numerically by engaging the MATLAB bvp4c function. The simultaneous influences of thermal radiation and Cattaneo–Christov heat flux (instead of conventional Fourier’s law of heat conduction) on the MHD nanofluid flow between two parallel plates were examined by Dogonchi and Ganji [21]. They found an analytical solution to the suggested problem by using the Duan–Rach approach (DRA), which facilitates finding undetermined coefficients without the help of numerical techniques. Sheikholeslami [22] deliberated the MHD non-Darcy nanofluid flow in a permeable cubic enclosure using the LBM (lattice Boltzmann method) and highlighted that the Nusselt number is an escalating function of permeability and buoyancy forces. Keeping in mind the applications of nanofluid flow in the attendance of magnetohydrodynamics, some featured recent explorations [23,24,25,26,27] have been added to the literature review.

The study of heat and mass transfer of nanofluid flows with chemical reactions past stretched surfaces has numerous significant applications in core industries such as chemical and metallurgical engineering. Examples include polymer production and food processing. Furthermore, heat and mass transfer models with chemical reaction effects are vital in several processes and have gained attention in recent years. Fascinating applications pertaining to chemical reactions embrace the distribution of temperature and moisture over agricultural fields and groves of fruit trees, the energy transfer of drying and evaporation at the surface of a water body, and the energy transfer of a wet cooling tower.

Entropy is pronounced the nonavailability of the system’s thermal energy for translation into mechanical work. The second law of thermodynamics asserts that throughout the process of conversion of energy into some beneficial work, there is forfeiture of energy that lowers the performance of energy conservation gadgets. Basically, entropy generation is directly proportionate to a loss of energy. Subsequently, entropy generation in a system results in the reduction of the quantity of energy present (exergy). Thus, the efficiency of the thermal system may be enhanced by slashing the entropy generation. In this regard, it is necessary to have the idea of energy generation’s distribution throughout a thermodynamic process with the intention of reduction in entropy generation. The idea of entropy generation was first floated by Bejan [28], who examined the reason for energy generation in a convective heat transfer model. The analysis of entropy generation in viscous fluid flow with an impact of suction and injection past a flat surface was conducted by Reveillere and Baytas [29]. The flow of MHD nanofluid past a spongy narrow vertical channel under the influences of nonlinear thermal radiation, entropy generation, and convective boundary conditions was considered by Lopez et al. [30]. Sheikholeslami [31] deliberated the analysis of exergy and entropy under the influence of Lorentz force past a permeable medium using a new numerical scheme, Control volume finite element method (CVFEM). Some useful recent explorations may be found in References [32,33,34,35,36].

The significance of non-Newtonian fluids cannot be denied due to their involvement in varied engineering and industrial applications such as petroleum production, the thinning of copper wires, and plastic manufacturing [37]. To exhibit such models, Naiver–Stokes equations are not enough since a single constitutive relation does not portray the core requirements of non-Newtonian fluids. That is why several discrete models have been proposed in the literature. Among these, viscoelastic fluid is the differential-type non-Newtonian fluid model that is simplest on its own and exhibits the characteristics of normal stress and shear rate viscosity. Hayat et al. [38] examined the flow of 3D viscoelastic nanofluid over an exponential stretched surface with mixed convection. The 3D viscoelastic nanofluid flow with nonlinear thermal radiation near a stagnation point was discussed by Farooq et al. [39]. Ramzan et al. [40,41,42] deliberated the flow of 3D viscoelastic nanofluid flow with the impacts of Newtonian heating, chemical reaction and MHD, and Soret–Dufour effects. The 3D flow of viscoelastic flow with temperature-dependent thermal conductivity was deliberated by Alsaedi et al. [43]. Some recent attempts considering viscoelastic fluids can be seen in References [44,45].

The literature review above indicates that few studies have been made on three-dimensional viscoelastic nanofluid flows associated with entropy generation and nonlinear thermal radiation past an exponential stretched surface. Therefore, we present our study on such an investigation by solving the proposed model analytically with three dimensionless distributions (Section 2) and examining the effects of different parameters (Section 5).

2. Mathematical Modeling

We considered the flow of MHD viscoelastic nanofluid past a permeable exponentially stretched surface with entropy generation. The analysis was performed in the attendance of nonlinear thermal radiation, chemical reaction, and heat generation and absorption with convective heat and mass boundary conditions. The magnetic field was applied parallel to the z axis. The Hall and electric field effects were overlooked. The induced magnetic field was also ignored due to our assumption of a small Reynolds number. The permeable surface was stretched in both the x and y directions with respective velocities $U_w$ and $V_w$, as shown in Figure 1.

The mathematical model, under assumptions, is represented by the following equations [46]:

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

(1)

$$\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}=v(\frac{{{\partial }^2}_u}{\partial z^2})-k_0(u\frac{{{\partial }^3}_u}{\partial x\partial z^2}+w\frac{{{\partial }^3}_u}{\partial z^3}-\frac{\partial u}{\partial x}\frac{{{\partial }^2}_u}{\partial z^2}-\frac{\partial u}{\partial z}\frac{{{\partial }^2}_w}{\partial z^2}-2\frac{\partial u}{\partial z}\frac{{{\partial }^2}_u}{dx\partial z}-2\frac{\partial w}{\partial z}\frac{{{\partial }^2}_u}{\partial z^2}]-\sigma \frac{\beta 0^2}{\rho }u,$$

(2)

$$u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}=v(\frac{{{\partial }^2}_v}{\partial z^2})-k_0(v\frac{{\partial }^{3}v}{\partial x\partial z^2}+w\frac{{{\partial }^3}_v}{\partial z^3}-\frac{\partial v}{\partial y}\frac{{{\partial }^2}_v}{\partial z^2}-\frac{\partial v}{\partial z}\frac{{{\partial }^2}_w}{\partial z^2}-2\frac{\partial v}{\partial z}\frac{{{\partial }^2}_v}{dy\partial z}-2\frac{\partial w}{\partial z}\frac{{{\partial }^2}_v}{\partial z^2})-\sigma \frac{\beta 0^2}{\rho }v,$$

(3)

$$u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z}=\alpha (\frac{{\partial }^{2}T}{\partial z^2})+\frac{Q}{\rho c_f}(T-T_{\infty })+\frac{(\rho c_p)}{(\rho c_f)}\left(D_B(\frac{\partial T}{\partial z}\frac{\partial C}{\partial z})+\frac{D_T}{T_{\infty }}{(\frac{\partial T}{\partial z})}^2\right)-\frac{1}{\rho c_p}\frac{\partial q_r}{\partial z},$$

(4)

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

(5)

subjected to the boundary conditions

$$\begin{array}{c}u=U_w={\mathrm{U}}_0{e}^{\frac{x+y}{L}}; v=V_w={\mathrm{V}}_0{e}^{\frac{x+y}{L}}, w=0,\\ -k\frac{\partial T}{\partial z}=h_f(T_f-T), -D_B\frac{\partial C}{\partial z}=h_s(C_s-C) at z=0\\ u\to 0, v\to 0, T\to T_{\infty }, C\to C_{\infty }\to \mathrm{as}\to z\to \infty .\end{array}$$

(6)

The surface stretching velocities, temperature, and concentration of the wall were given by

$$\begin{array}{c}{U}_w={\mathrm{U}}_0{e}^{\frac{x+y}{L}}; V_w={\mathrm{V}}_0{e}^{\frac{x+y}{L}}; T_w=T_f=T_{\infty }+T_0{e}^{\frac{A(x+y)}{L}}\\ C_w= C_s=C_{\infty }+C_0{e}^{\frac{B(x+y)}{L}}.\end{array}$$

(7)

The transformations are presented as given below:

$$\begin{array}{c}u = {\mathrm{U}}_0{e}^{\frac{x+y}{L}}{f}^{\prime }(\eta ), v={\mathrm{U}}_0{e}^{\frac{x+y}{L}}{g}^{\prime }(\eta ),w=-\sqrt{\frac{v{U}_0}{2L}}{e}^{\frac{x+y}{L}}(f+\eta f^{\prime }+g+\eta g^{\prime }),\\ T =T_{\infty }+{\mathrm{T}}_0{e}^{\frac{A(x+y)}{2L}} \theta (\eta ),C =C_{\infty }+{\mathrm{C}}_0{e}^{\frac{B(x+y)}{2L}} \varphi (\eta ),\eta =\sqrt{\frac{U_0}{2vL}}{e}^{\frac{x+y}{L}} z.\end{array}$$

(8)

Equation (1) was satisfied automatically, whereas Equations (2)–(5) took the form

$$f^{‴}+(f+g)f^{″}-2(f^{\prime }+g^{\prime })f^{\prime }-K(6f^{‴}{f}^{\prime }+(3g^{″}-3f^{″}+\eta g^{‴})f^{″}+(4g^{\prime }+2\eta g^{″})f^{‴})-(f+g+\eta g^{\prime })f^{⁗})-M^2{f}^{\prime }=0,$$

(9)

$$g^{‴}+(f+g)g^{″}-2(f^{\prime }+g^{\prime })g^{\prime }-K(6g^{‴}{g}^{\prime }+(3f^{″}-3g^{″}+\eta f^{‴})g^{″}+(4f^{\prime }+2\eta f^{″})g^{‴})-(f+g+\eta f^{\prime })f^{⁗}-M^2{g}^{\prime }=0,$$

(10)

$$\begin{array}{l}{\theta }^{″}+Pr((f+g){\theta }^{\prime }-A(f^{\prime }+g^{\prime })\theta +N_b{\theta }^{\prime }{\varphi }^{\prime }+N_t{{\theta }^{\prime }}^2+Rd(\frac{4}{3}{\theta }^{″}+4{({\theta }_w-1)}^2{\theta }^2{\theta }^{″}\\ +4({\theta }_w-1)\theta {\theta }^{″}+8({\theta }_w-1)\theta +8({\theta }_w-1)\theta {{\theta }^{\prime }}^2+4{{\theta }^{\prime }}^2+4{({\theta }_w-1)}^2{\theta }^2{{\theta }^{\prime }}^2\\ +\frac{4}{3}{({\theta }_w-1)}^3{\theta }^3{\theta }^{″}+PrQ\theta =0,\end{array}$$

(11)

$${\varphi }^{″}+LePr(f+g){\varphi }^{\prime }-BLePr(f^{\prime }+g^{\prime })\varphi +\frac{N_t}{N_b}{\theta }^{″}-PrLe\gamma \varphi =0,$$

(12)

subject to the transformed boundary conditions

$$\begin{array}{c}f=0, g=0, f^{\prime }=1,g^{\prime }=\beta , {\theta }^{\prime }=- {\lambda }_1(1-\theta (0)),{\varphi }^{\prime }=- {\lambda }_2(1-\varphi (0)),at \eta =0\\ f^{\prime }\to 0, g^{\prime }\to 0,\theta \to 0, \varphi \to 0 as \eta \to \infty .\end{array}$$

(13)

The parameters given in Equations (9)–(12) are defined as follows:

$$\begin{array}{c}K=\frac{k_0{U}_w}{2\nu L}, \beta =\frac{V_0}{U_0}, M^2=\frac{2\sigma B_0{}^{2}L}{\rho U_w}, Pr=\frac{\nu }{\alpha }, Nb=\frac{{(\rho c)}_p{D}_B(C_w-C_{\infty })}{{(\rho c)}_f\nu },\\ N_t=\frac{{(\rho c)}_p{D}_T(T_w-T_{\infty })}{{(\rho c)}_f\nu T_{\infty }},Le=\frac{\alpha }{D_B}, {\theta }_w=\frac{T_w}{T_{\infty }},Rd=\frac{4{\sigma }^{* }{T}_{\infty }{}^3}{3k^{*}{k}_1}, \gamma =\frac{h_f}{k}\sqrt{\frac{\nu }{\alpha }}.\end{array}$$

(14)

Skin Friction Coefficient and Local Nusselt and Sherwood Numbers

The equations of the skin friction coefficients along the x and y directions were given by

$$C_{fx}=\frac{{\tau }_{wx}}{1/2\rho U_w^2}=\frac{(\mu \frac{\partial u}{\partial z}+{\alpha }_1(u\frac{{{\partial }^2}_u}{\partial x\partial z}+v\frac{{{\partial }^2}_u}{\partial y\partial z}+w\frac{{{\partial }^2}_u}{\partial z^2}+\frac{\partial u}{\partial z}\frac{\partial u}{\partial x}+\frac{\partial v}{\partial z}\frac{\partial v}{\partial x}-\frac{\partial w}{\partial z}\frac{\partial u}{\partial z})){|}_{z=0} }{1/2\rho U_w^2},$$

(15)

$$C_{fy}=\frac{{\tau }_{wx}}{1/2\rho U_w^2}=\frac{(\mu \frac{\partial y}{\partial z}+{\alpha }_1(u\frac{{{\partial }^2}_v}{\partial x\partial z}+v\frac{{{\partial }^2}_v}{\partial y\partial z}+w\frac{{{\partial }^2}_v}{\partial z^2}+\frac{\partial u}{\partial z}\frac{\partial u}{\partial y}+\frac{\partial v}{\partial z}\frac{\partial v}{\partial y}-\frac{\partial w}{\partial z}\frac{\partial v}{\partial z})){|}_{z=0} }{1/2\rho U_w^2}.$$

(16)

The equations of the skin friction coefficients in dimensionless form were given by

$$C_{fx}={(\frac{Re}{2})}^{-\frac{1}{2}}(f^{″}+K(-(f+g)f^{‴}+5(f^{\prime }+g^{\prime })f^{″}+2f^{\prime }{f}^{″}+2g^{\prime }{g}^{″})){|}_{\eta =0},$$

(17)

$$C_{fy}={(\frac{Re}{2})}^{-\frac{1}{2}}{(}^{″}+K(-(f+g)f^{‴}+5(f^{\prime }+g^{\prime })g^{″}+2f^{\prime }{f}^{″}+2g^{\prime }{g}^{″})){|}_{\eta =0},$$

(18)

and the rates of heat and mass transfers in dimensionless forms were appended

$$N{u}_x=-\frac{x}{(T_w-T_{\infty })}(\frac{\partial T}{\partial z}){|}_{z=0}=-\frac{x}{L}{(\frac{Re}{2})}^{\frac{1}{2}}{\theta }^{\prime }(0),$$

(19)

$$S{h}_x=-\frac{x}{(C_w-C_{\infty })}(\frac{\partial C}{\partial z}){|}_{z=0}=-\frac{x}{L}{(\frac{Re}{2})}^{\frac{1}{2}}{\varphi }^{\prime }(0).$$

(20)

3. Convergence Analysis

The convergence of the series solutions was adjusted and regulated with the aid of auxiliary parameters ${\hslash }_{\mathrm{f}},{\hslash }_{\mathrm{g}},{\hslash }_{\mathsf{\theta }}$, and ${\hslash }_{\varnothing }$. For tolerable ranges, $\hslash$-curves were drawn and are given in Figure 2. These admissible ranges were ${\hslash }_{\mathrm{f}},{\hslash }_{\mathrm{g}},{\hslash }_{\mathsf{\theta }}$, and ${\hslash }_{\mathsf{\varphi }}$, which were $-1.85\ll {\hslash }_{\mathrm{f}}\ll -0.70,$ $-1.85\ll {\hslash }_{\mathrm{f}}\ll -0.70,$ $-1.80\ll {\hslash }_{\mathrm{g}}\ll -0.90,$ $-1.80\ll {\hslash }_{\mathsf{\theta }}\ll -0.70,$ and $-1.80\ll {\hslash }_{\mathsf{\varphi }}\ll -0.70.$ These values were corroborated by the numerical values given in

Table 1. Series solution convergence for varied orders of approximations when $\lambda$ = 0.1, Nb = 0.5, Nt = 0.3, K = 0.2, M = 0.3, Le = 1, Pr = 1.2, Rd = 0.2, A = B = 0.1.

Order of Approximation	$-{\mathit{f}}^{″}(0)$	$-{\mathit{g}}^{″}(0)$	$-{\mathit{\theta }}^{\prime }(0)$	$-\mathit{\varphi }(0)$
1	1.19588	0.12037	0.16113	0.16180
3	1.37115	0.13887	0.15382	0.15745
7	1.42903	0.14518	0.14911	0.15586
10	1.45762	0.14852	0.14206	0.15476
13	1.46000	0.14883	0.13937	0.15471
14	1.46000	0.14883	0.13798	0.15460
15	1.46000	0.14883	0.13796	0.15460
18	1.46000	0.14883	0.13796	0.15460

.

3.1. Homotopic Solutions

The homotopy analysis method necessitated the initial guesstimates ($f_0$, $g_0$, ${\theta }_0,{\varphi }_0$) with auxiliary linear operators ($L_f$, $L_g,L_{\theta },L_{\varphi }$) in the form [43]

$$f_0(\eta ) =1-e^{-\eta },g_0(\eta )=\mathsf{\lambda }(1-e^{-\eta }),{\theta }_0(\eta )=\frac{\gamma }{1+\gamma }{e}^{-\eta }, {\varphi }_0(\eta )=\frac{\gamma }{1+\gamma }{e}^{-\eta },$$

(21)

$$L_f=f^{‴}-f^{\prime },L_g=g^{‴}-g^{\prime },L_{\theta }={\theta }^{″}-\theta ,L_{\varphi }={\varphi }^{″}-\varphi$$

(22)

These auxiliary linear operators possessed the ensuing features

$$\begin{array}{c}{L}_f(B_1+B_2{e}^{\eta }+B_3{e}^{-\eta })=0,L_g(B_4+B_5{e}^{\eta }+B_6{e}^{-\eta })=0,\\ L_{\theta }(B_7{e}^{\eta }+B_8{e}^{-\eta })=0,L_{\varphi }(B_9{e}^{\eta }+B_{10}{e}^{-\eta })=0,\end{array}$$

(23)

where $B_i$, $i=1-10$ are arbitrary constants.

3.2. Deformation Problems at Zeroth Order

Here,

$$(1-p)L_f(\widehat{f}(\eta ;p)-f_0(\eta ))=p{\hslash }_f{N}_f(\widehat{f}(\eta ;p),\widehat{g}(\eta ;p)),$$

(24)

$$(1-p)L_g(\widehat{g}(\eta ;p)-g_0(\eta ))=p{\hslash }_g{N}_g(\widehat{f}(\eta ;p),\widehat{g}(\eta ;p)),$$

(25)

$$(1-p)L_{\theta }(\widehat{\theta }(\eta ;p)-{\theta }_0(\eta ))=p{\hslash }_{\theta }{N}_{\theta }(\widehat{f}(\eta ;p),\widehat{g}(\eta ;p),\widehat{\theta }(\eta ;p),\widehat{\varphi }(\eta ;p)),$$

(26)

$$(1-p)L_{\varphi }(\widehat{\varphi }(\eta ;p)-{\varphi }_0(\eta ))=p{\hslash }_{\mathsf{\varphi }}{N}_{\varphi }(\widehat{f}(\eta ;p),\widehat{g}(\eta ;p),\widehat{\theta }(\eta ;p),\widehat{\varphi }(\eta ;p)),$$

(27)

$$\widehat{g}(0;p)=0,{\widehat{f}}^{\prime }(\infty ;p)=0,{\widehat{g}}^{\prime }(\infty ;p)=0,$$

(28)

$$\widehat{g}(0;p)=0,{\widehat{g}}^{\prime }(0;p)=\beta ,{\widehat{g}}^{\prime }(\infty ;p)=0,$$

(29)

$$\widehat{\theta }(0;p)=-\gamma (1-\widehat{\theta }(0;p)),\widehat{\theta }(\infty ;p)=0,\varphi (0;p)=Nb\widehat{\varphi }(0;p)+Nt\widehat{\theta }(0;p),\varphi (\infty ;p)=0,$$

(30)

where ${\hslash }_{\mathrm{f}},{\hslash }_{\mathrm{g}},{\hslash }_{\mathsf{\theta }}$, and ${\hslash }_{\mathsf{\varphi }}$ characterize the nonzero convergence control parameters and $p\in [0,1]$ specifies the embedding parameter, whereas the nonlinear operators $N_f,N_g,N_{\theta }$, and $N_{\varphi }$ are specified by

$$\begin{gathered} N_f(\widehat{f}(\eta ;p),\widehat{g}(\eta ;p))=\frac{{\partial }^3\widehat{f}(\eta ;p)}{\partial {\eta }^3}-M^2\frac{\partial \widehat{f}(\eta ;p)}{\partial \eta }+(\widehat{f}(\eta ,p)+\widehat{g}(\eta ,p))\frac{{\partial }^2\widehat{f}(\eta ;p)}{\partial {\eta }^2}- \\ 2\left(\frac{\partial \widehat{f}(\eta ;p)}{\partial \eta }+\frac{\partial \widehat{g}(\eta ,p)}{\partial \eta }\right)\frac{\partial \widehat{f}(\eta ,p)}{\partial \eta }- \\ K(6\frac{\partial \widehat{f}(\eta ,p)}{\partial \eta }\frac{{\partial }^{3}f(\eta ,p)}{\partial {\eta }^3}+\left\{3\frac{{\partial }^2\widehat{g}(\eta ,p)}{\partial {\eta }^2}-3\frac{{\partial }^2\widehat{f}(\eta ,p)}{\partial {\eta }^2}-3\frac{{\partial }^2\widehat{f}(\eta ,p)}{\partial {\eta }^2}+\eta \frac{{\partial }^3\widehat{g}(\eta ,p)}{\partial {\eta }^3}\right\}\frac{{\partial }^2\widehat{f}(\eta ;p)}{\partial {\eta }^2}+ \\ \left\{4\frac{\partial \widehat{g}(\eta ,p)}{\partial \eta }+2\eta \frac{{\partial }^2\widehat{g}(\eta ,p)}{\partial {\eta }^2}\right\}\frac{{\partial }^3\widehat{f}(\eta ;p)}{\partial {\eta }^3})-\left(\widehat{f}(\eta ,p)+\widehat{g}(\eta ,p)+\eta \frac{\partial \widehat{g}(\eta ,p)}{\partial \eta }\right)\frac{{\partial }^4\widehat{f}(\eta ;p)}{\partial {\eta }^4}, \end{gathered}$$

(31)

$$\begin{gathered} N_g(\widehat{g}(\eta ;p),\widehat{g}(\eta ;p))=\frac{{\partial }^3\widehat{g}(\eta ;p)}{\partial {\eta }^3}-M^2\frac{\partial \widehat{g}(\eta ;p)}{\partial \eta }+(\widehat{f}(\eta ,p)+\widehat{g}(\eta ,p))\frac{{\partial }^2\widehat{g}(\eta ;p)}{\partial {\eta }^2}- \\ 2\left(\frac{\partial \widehat{f}(\eta ;p)}{\partial \eta }+\frac{\partial \widehat{g}(\eta ,p)}{\partial \eta }\right)\frac{\partial \widehat{g}(\eta ,p)}{\partial \eta }-K(6\frac{\partial \widehat{g}(\eta ,p)}{\partial \eta }\frac{{\partial }^3\widehat{g}(\eta ,p)}{\partial {\eta }^3}+ \\ \left\{3\frac{{\partial }^2\widehat{f}(\eta ,p)}{\partial {\eta }^2}-3\frac{{\partial }^2\widehat{g}(\eta ,p)}{\partial {\eta }^2}+\eta \frac{{\partial }^3\widehat{f}(\eta ,p)}{\partial {\eta }^3}\right\}\frac{{\partial }^2\widehat{g}(\eta ;p)}{\partial {\eta }^2}+4\left\{\frac{\partial \widehat{f}(\eta ,p)}{\partial \eta }+2\eta \frac{{\partial }^2\widehat{f}(\eta ,p)}{\partial {\eta }^2}\right\}\frac{{\partial }^3\widehat{g}(\eta ;p)}{\partial {\eta }^3})- \\ \left(\widehat{f}(\eta ,p)+\widehat{g}(\eta ,p)+\eta \frac{\partial \widehat{f}(\eta ,p)}{\partial \eta }\right)\frac{{\partial }^4\widehat{g}(\eta ;p)}{\partial {\eta }^4}, \end{gathered}$$

(32)

$$\begin{gathered} N_{\theta }(\widehat{f}(\eta ;p),\widehat{g}(\eta ;p),\widehat{\theta }(\eta ,p),\widehat{\varphi }(\eta ,p))=\frac{{\partial }^2\widehat{\theta }(\eta ,p)}{\partial {\eta }^2}+Pr(\widehat{f}(\eta ,p)+\widehat{g}(\eta ,p))\frac{\partial \widehat{\theta }(\eta ,p)}{\partial \eta } \\ -APr\left(\frac{\partial \widehat{f}(\eta ;p)}{\partial \eta }+\frac{\partial \widehat{g}(\eta ,p)}{\partial \eta }\right)\widehat{\theta }(\eta ,p)+NbPr\left(\frac{\partial \widehat{\theta }(\eta ,p)}{\partial \eta }\frac{\partial \widehat{\varphi }(\eta ,p)}{\partial \eta }\right)+NtPr{\left(\frac{\partial \widehat{\theta }(\eta ,p)}{\partial \eta }\right)}^2 \\ +\frac{4}{3}Rd{({\theta }_w-1)}^3{\widehat{\theta }}^3(\eta ,p)\frac{{\partial }^2\widehat{\theta }(\eta ,p)}{\partial {\eta }^2}+\frac{4}{3}Rd\frac{{\partial }^2\widehat{\theta }(\eta ,p)}{\partial {\eta }^2}+4Rd({\theta }_w-1){\widehat{\theta }}^2(\eta ,p)\frac{{\partial }^2\widehat{\theta }(\eta ,p)}{\partial {\eta }^2} \\ +4Rd({\theta }_w-1)\widehat{\theta }(\eta ,p)\frac{{\partial }^2\widehat{\theta }(\eta ,p)}{\partial {\eta }^2}+8Rd({\theta }_w-1)\widehat{\theta }(\eta ,p){\left(\frac{\partial \widehat{\theta }(\eta ,p)}{\partial \eta }\right)}^2+ \\ +4Rd{\left(\frac{\partial \widehat{\theta }(\eta ,p)}{\partial \eta }\right)}^2+4Rd{({\theta }_w-1)}^2{\widehat{\theta }}^2(\eta ,p){\left(\frac{\partial \widehat{\theta }(\eta ,p)}{\partial \eta }\right)}^2+PrQ\widehat{\theta }(\eta ,p), \end{gathered}$$

(33)

$$\begin{gathered} N_{\varphi }(\widehat{f}(\eta ;p),\widehat{g}(\eta ;p),\widehat{\theta }(\eta ,p),\widehat{\varphi }(\eta ,p))=\frac{{\partial }^2\widehat{\varphi }(\eta ,p)}{\partial {\eta }^2}+LePr\left(\widehat{f}(\eta ,p)+\widehat{g}(\eta ,p)\right) \\ \frac{\partial \widehat{\varphi }(\eta ,p)}{\partial \eta }-BLePr(\frac{\partial \widehat{f}(\eta ,p)}{\partial \eta }+\frac{\partial \widehat{g}(\eta ,p)}{\partial \eta })\widehat{\varphi }(\eta ,p)+\frac{Nt}{Nb}\frac{{\partial }^2\widehat{\theta }(\eta ,p)}{\partial {\eta }^2} \\ +\frac{Nt}{Nb}\frac{{\partial }^2\widehat{\theta }(\eta ,p)}{\partial {\eta }^2}-PrLe\widehat{\varphi }(\eta ,p). \end{gathered}$$

(34)

For $p=0$ and $p=1$, we get

$$\begin{array}{l}\widehat{f}(\eta ;0)=f_0(\eta ),\widehat{g}(\eta ;0)=g_0(\eta ),\widehat{\theta }(\eta ;0)={\theta }_0(\eta ),\widehat{\varphi }(\eta ;0)={\varphi }_0(\eta )-\gamma (1-\widehat{\theta }(\eta ;0)),\\ \widehat{\varphi }(\eta ;0)=Nb\widehat{\varphi }(\eta ;0)+Nt\widehat{\theta }(\eta ;0),\widehat{f}(\eta ;1)=f(\eta ),\widehat{g}(\eta ;1)=g(\eta ),\widehat{\theta }(\eta ;1)=\theta (\eta ),\widehat{\varphi }(\eta ;1)=\varphi (\eta )\end{array}$$

(35)

The subsequent expressions were attained via Taylor’s series expansion:

$$\widehat{f}(\eta ;p)=f_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }{f}_m{p}^m},f_m(\eta )={\frac{1}{m!}\frac{{\partial }^m\widehat{f}(\eta ;p)}{\partial p^m}|}_{p=0},$$

(36)

$$\widehat{g}(\eta ;p)=g_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }{g}_m{p}^m},g_m(\eta )={\frac{1}{m!}\frac{{\partial }^m\widehat{g}(\eta ;p)}{\partial p^m}|}_{p=0},$$

(37)

$$\widehat{\theta }(\eta ;p)={\theta }_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }{\theta }_m}(\eta )p^m,{\theta }_m(\eta )={\frac{1}{m!}\frac{{\partial }^m\widehat{\theta }(\eta ;p)}{\partial p^m}|}_{p=0},$$

(38)

$$\widehat{\varphi }(\eta ;p)={\varphi }_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }{\varphi }_m}(\eta )p^m,{\varphi }_m(\eta )={\frac{1}{m!}\frac{{\partial }^m\varphi (\eta ;p)}{\partial p^m}|}_{p=0}.$$

(39)

The convergence of control parameters ${\hslash }_{\mathrm{f}},{\hslash }_{\mathrm{g}}, {\hslash }_{\mathsf{\theta }}$, and ${\hslash }_{\mathsf{\varphi }}$ were selected in such a way that the series (36)–(39) converged at $p=1$. Then

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

(40)

$$\widehat{g}(\eta ;p)=g_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }{g}_m}(\eta ),$$

(41)

$$\widehat{\theta }(\eta ;p)={\theta }_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }{\theta }_m}(\eta ),$$

(42)

$$\widehat{\varphi }(\eta ;p)={\varphi }_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }{\varphi }_m}(\eta ).$$

(43)

3.3. The m-th Order Problem

The equations of the m-th order are

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

(44)

$$L_g(g_m(\eta )-X_m{g}_{m-1}(\eta ))={\hslash }_g{R}_g^m(\eta ),$$

(45)

$$L_{\theta }({\theta }_m(\eta )-X_m{\theta }_{m-1}(\eta ))={\hslash }_{\theta }{R}_{\theta }^m(\eta ),$$

(46)

$$L_{\varphi }({\varphi }_m(\eta )-X_m{\varphi }_{m-1}(\eta ))={\hslash }_{\mathsf{\varphi }}{R}_{\varphi }^m(\eta ),$$

(47)

$$\begin{gathered} R_f^m(\eta )=f_{m-1}^{‴}(\eta )-M^2{f}_{m-1}^{\prime }(\eta )+{\displaystyle \sum _{k=0}^{m-1}(f_{m-1-k}{f}_k^{″}+g_{m-1-k}^{\prime }{f}_k^{″})}-2{\displaystyle \sum _{k=0}^{m-1}{f}_{m-1-k}^{\prime }}{f}_k^{\prime } \\ -2{\displaystyle \sum _{k=0}^{m-1}{g}_{m-1-k}^{\prime }}{f}_k^{\prime }-k(6{\displaystyle \sum _{k=0}^{m-1}{f}_{m-1-k}^{\prime }}{f}_k^{‴}+3{\displaystyle \sum _{k=0}^{m-1}{g}_{m-1-k}^{″}}{f}_k^{″}-3{\displaystyle \sum _{k=0}^{m-1}{f}_{m-1-k}^{″}}{f}_k^{″}+{\displaystyle \sum _{k=0}^{m-1}\eta g_{m-1-k}^{‴}}{f}_k^{″} \\ +4{\displaystyle \sum _{k=0}^{m-1}{g}_{m-1-k}^{\prime }}{f}_k^{‴}+2{\displaystyle \sum _{k=0}^{m-1}\eta g_{m-1-k}^{″}}{f}_k^{‴}-{\displaystyle \sum _{k=0}^{m-1}{f}_{m-1-k}}{f}_k^{⁗}-{\displaystyle \sum _{k=0}^{m-1}\eta g_{m-1-k}^{\prime }}{f}^{⁗}), \end{gathered}$$

(48)

$$\begin{gathered} R_g^m(\eta )=g_{m-1}^{‴}(\eta )-M^2{g}_{m-1}^{\prime }(\eta )+{\displaystyle \sum _{k=0}^{m-1}(f_{m-1-k}{g}_k^{″}+g_{m-1-k}^{\prime }{g}_k^{″})}-2{\displaystyle \sum _{k=0}^{m-1}{g}_{m-1-k}^{\prime }}{g}_k^{\prime } \\ -2{\displaystyle \sum _{k=0}^{m-1}{g}_{m-1-k}^{\prime }}{f}_k^{\prime }-k(6{\displaystyle \sum _{k=0}^{m-1}{g}_{m-1-k}^{\prime }}{f}_k^{‴}+3{\displaystyle \sum _{k=0}^{m-1}{f}_{m-1-k}^{″}}{g}_k^{″}-3{\displaystyle \sum _{k=0}^{m-1}{g}_{m-1-k}^{″}}{g}_k^{″}+{\displaystyle \sum _{k=0}^{m-1}\eta f_{m-1-k}^{‴}}{g}_k^{″} \\ +4{\displaystyle \sum _{k=0}^{m-1}{f}_{m-1-k}^{\prime }}{g}_k^{‴}+2{\displaystyle \sum _{k=0}^{m-1}\eta f_{m-1-k}^{″}}{g}_k^{‴}-{\displaystyle \sum _{k=0}^{m-1}{g}_{m-1-k}}{g}_k^{⁗}-{\displaystyle \sum _{k=0}^{m-1}\eta f_{m-1-k}^{\prime }}{g}^{⁗}), \end{gathered}$$

(49)

$$\begin{gathered} R_{\theta }^m(\eta )={\theta }_{m-1}^{″}(\eta )+Pr{\displaystyle \sum _{k=0}^{m-1}(f_{m-1-k}{\theta }_k^{\prime }+g_{m-1-k}{\theta }_k^{\prime })}-APr{\displaystyle \sum _{k=0}^{m-1}(f_{m-1-k}^{\prime }{\theta }_k+g_{m-1-k}^{\prime }{\theta }_k)} \\ +NbPr{\displaystyle \sum _{k=0}^{m-1}{\theta }_{m-1-k}^{\prime }}{\varphi }_k^{\prime }+NtPr{\displaystyle \sum _{k=0}^{m-1}{\theta }_{m-1-k}^{\prime }}{\theta }_k^{\prime }+Rd(\frac{4}{3}{({\theta }_w-1)}^3{\displaystyle \sum _{k=0}^{m-1}{\theta }_{m-1-k}^3}{\theta }_k^{″}+ \\ \frac{4}{3}{\theta }_{m-1}^{″}(\eta )+4{({\theta }_w-1)}^2{\displaystyle \sum _{k=0}^{m-1}{\theta }_{m-1-k}^2}{\theta }_k^{″}+4({\theta }_w-1){\displaystyle \sum _{k=0}^{m-1}{\theta }_{m-1-k}}{\theta }_k^{″}+ \\ 8({\theta }_w-1){\displaystyle \sum _{k=0}^{m-1}{\theta }_{m-1-k}}{\theta }_k^{″}+4{\theta }_{m-1}^{\prime }(\eta )+4{({\theta }_w-1)}^2{\displaystyle \sum _{k=0}^{m-1}{\theta }_{m-1-k}^2}{\theta }_k^{{\prime }^2})+PrQ{\theta }_{m-1}(\eta ), \end{gathered}$$

(50)

$$R_{\varphi }^m(\eta )={\varphi }_{m-1}^{″}(\eta )+LePr{\displaystyle \sum _{k=0}^{m-1}(f_{m-1-k}{\varphi }_k^{\prime }-g_{m-1-k}{\varphi }_k^{\prime })}-BLePr,$$

(51)

$$X_m=\{\begin{array}{l}0m\le 1\\ 1m>1\end{array}.$$

(52)

The general solutions $(f_m,g_m,{\theta }_m,{\varphi }_m)$ for Equations (48)–(51) in terms of special solutions $({f^{*}}_m,{g^{*}}_m,{{\theta }^{*}}_m,{{\varphi }^{*}}_m)$ are given by

$$f_m(\eta )=f_m^{*}(\eta )+C_1+C_2{e}^{\eta }+C_3{e}^{-\eta },$$

(53)

$$g_m(\eta )=g_m^{*}(\eta )+C_4+C_5{e}^{\eta }+C_6{e}^{-\eta },$$

(54)

$${\theta }_m(\eta )={\theta }_m^{*}(\eta )+C_7{e}^{\eta }+C_8{e}^{-\eta },$$

(55)

$${\varphi }_m(\eta )={\varphi }_m^{*}(\eta )+C_9{e}^{\eta }+C_{10}{e}^{-\eta }.$$

(56)

4. Entropy Analysis

The equation of entropy generation [47,48,49] was given by

$$\begin{array}{ll}{S}_G& =\underset{a}{\underbrace{\frac{K}{{T_{\infty }}^2}\left(1+\frac{16{\sigma }^{*}{T}^3}{3K{K}^{*}}\right){(\frac{\partial T}{\partial z})}^2}}+\underset{b}{\underbrace{\frac{R{D}_B}{C_{\infty }}{(\frac{\partial C}{\partial z})}^2+\frac{R{D}_B}{T_{\infty }}(\frac{\partial C}{\partial z}\frac{\partial T}{\partial z})}}\\ & +\underset{c}{\underbrace{\frac{\rho C_p}{T_{\infty }}D}}+\underset{d}{\underbrace{\frac{\sigma {B_0}^2}{T_{\infty }}(u^2+v^2)}},\end{array}$$

(57)

where a, b, c and d represent the thermal irreversibility or entropy generation due to heat transfer, the concentration irreversibility or entropy generation due to mass transfer, the entropy generation due to viscous effects in the fluid, and the entropy generation due to magnetic effects of the fluid, respectively. Then

$$D={(\frac{\partial C}{\partial z})}^2\left\{\begin{array}{l}(\frac{\partial u}{\partial x})(\frac{2\nu }{C_p(1+{\lambda }_1)}(\frac{\partial u}{\partial x})+\frac{2\nu {\lambda }_2}{C_p(1+{\lambda }_1)}(u\frac{{\partial }^{2}u}{\partial x^2}+v\frac{{\partial }^{2}u}{\partial x\partial y}+w\frac{{\partial }^{2}u}{\partial x\partial z}))+\\ (\frac{\partial v}{\partial y})(\frac{2\nu }{C_p(1+{\lambda }_1)}(\frac{\partial v}{\partial x})+\frac{2\nu {\lambda }_2}{C_p(1+{\lambda }_1)}(u\frac{{\partial }^{2}v}{\partial x\partial y}+v\frac{{\partial }^{2}v}{{\partial }^{2}y}+w\frac{{\partial }^{2}v}{\partial y\partial z}))\end{array}\right\},$$

(58)

and

$${(S_{gen})}_0=\frac{K}{T_{\infty }{}^2}{\left(\frac{∆T}{L^2}\right)}^2$$

(59)

is the characteristic entropy generation rate. Using the transformations in Equation (6), Equation (21) takes the form

$$\begin{array}{ll}Ns=& \frac{S_{gen}}{{(S_{gen})}_0}=\frac{Re}{2}{e}^{\chi }\left(\left(1+\frac{4Rd}{3}{(\theta ({\theta }_w-1)+1)}^3\right){{\theta }^{\prime }}^2+\frac{{\sum }^2\epsilon }{{\mathsf{\Omega }}^2}{{\varphi }^{\prime }}^2+\frac{\sum \epsilon }{\mathsf{\Omega }}{\theta }^{\prime }{\varphi }^{\prime }\right)\\ & +\frac{Br{e}^{\chi }}{\mathsf{\Omega }(1+{\lambda }_1)}d+\frac{MBr}{\mathsf{\Omega }}{f^{\prime }}^2\end{array}$$

(60)

where

$$\begin{array}{l}d=(2f^{\prime }+\eta f^{″})\left((2f^{\prime }+\eta f^{″})+\frac{\beta S}{2}((4f^{\prime }+2\eta f^{″})(f^{\prime }+g^{\prime })-(3f^{″}+\eta f^{‴})(f+g))\right)+\\ (2g^{\prime }+\eta g^{″})((2g^{\prime }+\eta g^{″})+\frac{\beta S}{2}((4g^{\prime }+2\eta g^{″})(f^{\prime }+g^{\prime })-(3g^{″}+\eta g^{‴})(f+g)))\end{array}$$

(61)

An important Bejan number defined by the ratio of entropy generation due to thermal effect and total entropy generated by thermal, concentration, and fluid frictions forces was as follows [50]:

$$Be=\frac{\mathrm{The} \mathrm{entropy} \mathrm{generation} \mathrm{due} \mathrm{to} \mathrm{thermal} \mathrm{irreversibility}}{\mathrm{The} \mathrm{total} \mathrm{entropy} \mathrm{generation}}.$$

(62)

In dimensionless form, it became

$$Be=\frac{\frac{Re}{2}{e}^{\chi }\left(1+\frac{4Rd}{3}{(\theta ({\theta }_w-1)+1)}^3\right){{\theta }^{\prime }}^2}{\frac{Re}{2}{e}^{\chi }\left(\left(1+\frac{4Rd}{3}{(\theta ({\theta }_w-1)+1)}^3\right){{\theta }^{\prime }}^2+\frac{{\sum }^2\epsilon }{{\mathsf{\Omega }}^2}{{\varphi }^{\prime }}^2+\frac{\sum \epsilon }{\mathsf{\Omega }}{\theta }^{\prime }{\varphi }^{\prime }\right)+\frac{Br{e}^{\chi }}{\mathsf{\Omega }(1+{\lambda }_1)}d+\frac{MBr}{\mathsf{\Omega }}{f^{\prime }}^2}.$$

(63)

The value of the irreversibility parameter must lie between $0<Be<1$. if $Be=0,$ then there is no entropy generated due to heat transfer. When $Be<0.5$, then the entropy generation due to fluid friction is dominant over the entropy generation due to heat transfer. On the other hand, if $Be>0$, then reverse relation is obvious for the entropy generation due to heat and total entropy generation.

The averaged entropy generation number could be evaluated using the following integral formula:

$${[Ns]}_{avg}=\frac{1}{\forall }{\displaystyle \underset{\forall }{\int }Nsd\forall },$$

(64)

where $\forall$ is the length of the boundary layer region.

5. Results and Discussion

This section is dedicated to having a deep insight into changes of sundry parameters on involved profiles. The ranges of all parameters were selected via Turkyilmazoglu [51]. Figure 3 was drawn to show the impact of viscoelastic parameter $K$ on both velocities. Both velocity distributions were declining functions of $K$. Augmented values of the $K$ meant the development of tensile stress, which ultimately lowered both velocity components. The impression of magnetic parameter $M$ on both velocity components is depicted in Figure 4. A diminishing behavior was witnessed for both velocity components. A strong Lorentz force was perceived due to the magnetic field, hindering the movement of the fluid flow. Thus, the reduction in velocity components was noticed. Figure 5 is outlined to grasp the impact of chemical reaction parameter $\gamma$ and Prandtl number $Pr$ on the concentration field. It was visualized that the rate of mass transfer was on the decline with augmented values of $\mathsf{\gamma }.$ Thus, the respected boundary layer was also enhanced. Smaller values of thermal diffusivity augmented the $Pr$ and thereby lowered the concentration distributions. The outcome of temperature exponent $A$ and Prandtl number $Pr$ on temperature profile can be visualized in Figure 6. The Prandtl number is the quotient of momentum to thermal diffusivity. Smaller numbers for thermal diffusivity meant higher values of $Pr$, thus lowering the temperature. A similar effect was noted in the case of $A$. Figure 7 was illustrated to depict the influence of thermophoretic parameter $Nt$ and Brownian motion parameter $Nb$ on temperature field. It was comprehended that temperature was an increasing function of both parameters. This is because smaller particles were pushed toward the colder surface from the hotter one. Another reason for this phenomenon is that enhanced random movement of the fluid’s particles also became the main source of heightening the temperature of the fluid. To understand the influence of Brownian motion parameter $Nb$ and thermophoresis parameter $Nt$ on the concentration distribution, Figure 8 was plotted. Higher values of both parameters augmented the concentration profile. Smaller particles were pulled toward the colder region due to increased values of the thermophoretic parameter, thus lowering the concentration of the fluid. Figure 9 was drawn to grasp the impact of radiation parameter $Rd$ and heat absorption parameter $Q$ on the temperature field. In the experiment, the temperature rose with higher values of $Rd$. This is because growing values of Rd meant more emission of heat that ultimately raised the temperature of the fluid. Also, the higher estimates of $Q$ produced more heat that also resulted in an increase in the fluid’s temperature. The variation of magnetic parameter M for the entropy generation $Ns$ is presented in Figure 10. Initially, due to drag forces produced by the magnetic field, the energy was produced, which increased the entropy of the fluid. However, as it moved away from the plate, the influence of these forces reduced gradually and eventually lessened the effect of entropy as well. The effect of viscoelastic parameter $K$ on entropy generation and the Bejan number $Be$ is displayed in Figure 11 and Figure 12, respectively. Both the figures clearly indicate that the enhancement of values of $K$ enhanced entropy generation as well as the Bejan number near the surface due to an increase in viscous effects. However, the opposite trend was noticed in both cases while moving away from the surface. In Figure 13, a relationship between the thermal radiation parameter $Rd$ and dimensionless Bejan number $Be$ is portrayed. An increment in thermal radiation parameter $Rd$ enhanced the internal heat generation in the moving fluid, which caused the Bejan number to escalate. Figure 14, Figure 15, Figure 16 and Figure 17 were illustrated to show the impacts of dimensionless temperature difference Ω, Reynolds number Re, thermal radiation parameter Rd, and Brinkman number Br on the variation of local entropy generation number Ns, respectively. From Figure 14, it is revealed that as the temperature difference increased, a decrease in the entropy generation parameter was noticed. Larger values of Ω did not affect the viscous part of the entropy generation. However, only the thermal entropy generation segment was responsible for higher temperature gradients in the transverse direction. In Figure 15, Figure 16 and Figure 17, Ns was an increasing function of these three: Re, Rd, and Br. Figure 18 and Figure 19 were drawn to depict the influence of M, K, Pr, and Rd on average entropy Ns_(avg). From these figures, it was gathered that Ns_(avg) was an escalating function of all four parameters. This meant irreversibility mounted as the group parameters grew.

Table 2. Values of local Nusselt and Sherwood numbers for involved parameters.

λ	Nb	Nt	Le	Pr	M	K	Rd	A	Q	$-{(\frac{\mathit{R}\mathit{e}}{2})}^{-\frac{1}{2}}\mathit{N}{\mathit{u}}_{\mathit{x}}.$	$-{(\frac{\mathit{R}\mathit{e}}{2})}^{-\frac{1}{2}}\mathit{S}{\mathit{h}}_{\mathit{x}}$
0.1	-	-	-	-	-	-	-	-	-	0.13878	0.12199
0.2	-	-	-	-	-	-	-	-	-	0.1388	0.12157
0.5	-	-	-	-	-	-	-	-	-	0.1421	0.12032
-	0.5	-	-	-	-	-	-	-	-	0.13878	0.12199
-	1.0	-	-	-	-	-	-	-	-	0.13772	0.11917
-	1.5	-	-	-	-	-	-	-	-	0.13770	0.11823
-	-	0.0	-	-	-	-	-	-	-	0.13878	0.11635
-	-	0.2	-	-	-	-	-	-	-	0.13878	0.12199
-	-	0.5	-	-	-	-	-	-	-	0.13878	0.12576
-	-	-	1.0	-	-	-	-	-	-	0.14447	0.12981
-	-	-	1.5	-	-	-	-	-	-	0.13878	0.12199
-	-	-	2.0	-	-	-	-	-	-	0.12572	0.11010
-	-	-	-	1.0	-	-	-	-	-	0.13878	0.12199
-	-	-	-	1.2	-	-	-	-	-	0.13775	0.08889
-	-	-	-	1.5	-	-	-	-	-	0.13774	0.06499
-	-	-	-	-	0.0	-	-	-	-	0.13878	0.12199
-	-	-	-	-	0.2	-	-	-	-	0.13878	0.12199
-	-	-	-	-	0.3	-	-	-	-	0.13878	0.12199
-	-	-	-	-	-	0.0	-	-	-	0.13878	0.12199
-	-	-	-	-	-	0.02	-	-	-	0.13878	0.12199
-	-	-	-	-	-	0.04	-	-	-	0.13878	0.12199
-	-	-	-	-	-	-	0.2	-	-	0.13878	0.12199
-	-	-	-	-	-	-	0.4	-	-	0.14392	0.12146
-	-	-	-	-	-	-	0.5	-	-	0.14480	0.12123
-	-	-	-	-	-	-	-	0.1	-	0.13878	0.12199
-	-	-	-	-	-	-	-	0.5	-	0.13878	0.12199
-	-	-	-	-	-	-	-	0.7	-	0.13878	0.12199
-	-	-	-	-	-	-	-	-	0.2	0.13878	0.12199
-	-	-	-	-	-	-	-	-	0.4	0.14962	0.12199
-	-	-	-	-	-	-	-	-	0.5	0.15679	0.12199

was initiated for the rate of heat and mass transfers versus varied involved parameters. It was noticed that the Nusselt number was a growing function of $\lambda ,Rd, Q$ and was a declining function of $Nb, Le, Pr.$ Similarly, the Sherwood number decreased for the values of $\lambda , Nb, Le, Pr, Rd,$ and increased for the values of $Nt.$

Table 3. Values of skin friction coefficients for involved parameters.

λ	M	K	$-{(\frac{\mathit{R}\mathit{e}}{2})}^{1/2}{\mathit{C}}_{{\mathit{f}}_{\mathit{x}}}$	$-{(\frac{\mathit{R}\mathit{e}}{2})}^{1/2}{\mathit{C}}_{{\mathit{f}}_{\mathit{y}}}$
0.1	-	-	1.6768	0.2237
0.2	-	-	1.7698	0.4089
0.5	-	-	2.0571	1.0804
-	0.3	-	1.6768	0.2237
-	0.5	-	1.7422	0.2325
-	1.0	-	2.0212	0.2697
-	-	0.02	1.6768	0.2237
-	-	0.03	1.9607	0.2675
-	-	0.04	1.8168	0.3138

depicts the values of both skin friction coefficients along $\mathrm{the} x$ and $y$ directions. It was gathered that both coefficients were growing functions of $\lambda , M,$ and $K.$

Table 4. Comparison of present values to Liu et al. [52] in the limiting case when $K=Nb=Nt=Rd=Q=M=0$ (also, values of the convective boundary were neglected).

β	Pr	A	Liu et al. [52]	Present Study
0.0	0.7	0.0	−0.42583804	−0.4258120
		2.0	−1.02143617	−1.0214514
		5.0	−1.64165922	−1.6416620
0.25	0.7	0.0	−0.47609996	−0.4761032
		2.0	−1.14199997	−1.1420014
		5.0	−1.83543073	−1.8354210
0.50	0.7	0.0	−0.52154103	−0.5215267
		2.0	−1.25099820	−1.2509991
		5.0	−2.01061361	−2.0106021
0.75	0.7	0.0	−0.56332861	−0.5633148
		2.0	−1.35123246	−1.3512221
		5.0	−2.17171091	−2.1717006
1.0	0.7	0.0	−0.60222359	−0.6022167
		2.0	−1.44452826	−1.4445214
		5.0	−2.32165661	−2.3216340

was erected to check the reliability of the presented solution by comparing it to Liu et al. [52] in the limiting case. An excellent concurrence was obtained, confirming the reliability of the presented results.

6. Concluding Remarks

Entropy generation of 3D viscoelastic nanofluid fluid flow past an exponential stretched surface with nonlinear thermal radiation, chemical reaction, and magnetohydrodynamics was discussed analytically with the homotopy analysis method. The whole analysis was completed with the effect of heat generation and absorption supported by convective heat and mass boundary conditions. The main findings of the analysis are:

• The velocity components were declining functions of the viscoelastic parameter.
• The temperature field improved with an increase in radiation parameter.
• Thermophoresis and Brownian motion parameters had an opposite effect on concentration distribution.
• With growing values of the magnetic parameter, both velocity components declined.
• The Bejan number is an increasing function of the thermal radiation parameter.
• Entropy generation decreased for escalating values of the temperature difference parameter.