Numerical Investigation of Natural Convection Viscoelastic Jeffrey’s Nanofluid Flow from a Vertical Permeable Flat Plate with Heat Generation, Thermal Radiation, and Chemical Reaction

The boundary layer flow of an incompressible viscoelastic Jeffrey’s nanofluid from a vertical permeable flat plate is investigated. We consider the effects of heat generation, thermal radiation, and chemical reaction on the fluid flow. The nonlinear transformed coupled differential equations that describe the transport processes are solved numerically using a multidomain bivariate spectral quasilinearization method (MD-BSQLM). This innovative method involves blending the quasilinearization idea with the bivariate Lagrange interpolation. The solutions of the resulting system of equations are then obtained sequentially on multiple intervals using the Chebyshev spectral collocation method. The method is shown to give accurate solutions for boundary layer-type equations. The influence of various physical parameters on velocity, temperature, and nanoparticle concentration fields, as well as on the skin friction and heat and mass transfer coefficients, is shown and discussed in detail. The range of the values of the governing parameters considered in this study is between 1⁄20, 4 . For qualitative validation of the results and the numerical method used, calculations were carried out to graphically obtain the velocity, temperature, and nanoparticle concentration fields for selected physical parameter values. The results obtained were found to correlate with the results from published literature. For quantitative verification of our findings, the MD-BSQLM numerical solutions were again confirmed against published results reported in the literature, and the results were observed to be in perfect agreement. This study’s findings indicate that the Deborah number and suction parameter have related effects on the velocity profile, which is to suppress both the flow velocity and the momentum boundary layer thickness. Increasing the heat generation and thermal radiation parameters enhances both the temperature and thermal boundary layer depths. In contrast, an increase in the chemical reaction parameter causes a decrease in the fluid concentration.


Introduction
The study of non-Newtonian fluid flow has been an active area of interest in the past several decades, specifically concerning heat and mass transfer processes in industrial processes such as the manufacture of plastic films and artificial fibers, cooling of metallic sheets, aerodynamic extrusion of plastic sheets, liquid film condensation, and crystal growing. According to Shehzad et al. [1], non-Newtonian fluids do not obey Newton's law of viscosity. In such fluids, nonlinearity exists between the shear stress and the strain rate relation, thereby making the flow model complicated and relating the shear stresses to the velocity field [2]. Significant contributions to the study of different non-Newtonian fluid models by different researchers include oblique flows [3], nanofluid flow [4], Maxwell fluid [5], tangent hyperbolic flows [6], Williamson fluid [7], Khan et al. [8], Bibi et al. [9], Powell-Eyring fluid [10], Oldroyd-B fluid [11], power-law fluid [12], Walter's B fluid [13], and Jeffrey's fluid [14]. An interesting and one of the simplest subclasses of non-Newtonian fluids is Jeffrey's fluid. This is because Jeffrey's fluid model utilizes the time derivatives instead of converted derivatives and degenerates to the Newtonian model at very high wall shear stress [2]. Also, Jeffrey's fluid model approximates well the rheological behavior of a wide range of industrial processes such as biotechnological detergents, physiological suspensions, dense foams, geological sediments, and cosmetic creams. Due to the importance of Jeffrey's fluids, many researchers have studied several flow geometries of this fluid with different boundary conditions, using various techniques to solve such models. These include the study of Turkyilmazoglu and Pop [15] who presented an analytical solution for the flow and heat transfer of Jeffrey's fluid near the stagnation point on a stretching/shrinking sheet with the parallel external flow. In their study, it was shown that the structure of the solutions they presented strongly depends on the stretching strength parameter; this parameter measures the ratio of the strength of the external flow to surface stretching/shrinking. When this parameter is set to zero, the solutions evolve into the multiple/triple solutions already known in the literature. The influence of heat transfer on peristaltic transport of Jeffrey's fluid in a vertical porous stratum was studied by Vajravelu et al. [16] using the perturbation technique. In their study, it was observed that the effects of Jeffrey's number, the Grashof number, the perturbation parameter, and the peristaltic wall deformation parameter are the strongest on the trapping bolus phenomenon.
The combined effect of heat and mass transfer on Jeffrey's fluid over a stretching sheet in the presence of a heat source/heat sink was investigated by Qasim [17]. His study assumed that the surface temperature and the concentration vary according to the power-law form. Exact solutions for the nonlinear equations governing the flow were derived by the power series method using Kummer's confluent hypergeometric functions. He further observed that the velocity increases with an increase in the Deborah number. Simultaneously, the temperature is a decreasing function of the Deborah number, and the thermal boundary layer thickness decreases by increasing the wall temperature and heat sink parameters. Ahmad et al. [18] numerically investigated the steady two-dimensional magnetohydrodynamic (MHD) mixed convection boundary layer flow and heat transfer of Jeffrey's fluid over an exponentially stretched plate using an implicit finite difference scheme. In their investigation, local similarity solutions were obtained for some embedded parameters, such as the Deborah number, mixed convection parameter, Prandtl number, and Hartmann number. They also observed that as the Deborah number increases, unusual behavior is found to happen when the flow is assisted by the buoyancy force. The shooting method with the fourth-order Runge-Kutta scheme was used by Narayana and Babu [19] to study the effects of chemical reaction and a heat source on MHD heat and mass transfer of an electrically conducting Jeffrey's fluid over a stretching sheet in the presence of a power-law form of temperature and concentration. It was observed in their study that the Deborah number and ratio of relaxation and retardation time parameter have opposite effects on the skin friction coefficient. Saqib et al. [20] reported the applications of Caputo-Fabrizio timefractional derivatives to generalize Jeffrey's fluid over a vertical static plate. Analytical solutions using the Laplace transform technique was presented for the governing equations of generalized Jeffrey's fluid model. They also reported that free convection is caused due to the temperature gradient; therefore, heat transfer was considered for free convection in their study. Exact solutions for unsteady free convection flow of Jeffrey's fluid were obtained by Khan [21] using the Laplace transform technique. His research found that the obtained solutions satisfy the imposed initial and boundary conditions and can be easily reduced to similar solutions for a Newtonian fluid. Zin et al. [22] presented exact and numerical solutions for unsteady heat and mass transfer problems of Jeffrey's fluid with MHD and Newtonian heating effects using the Laplace transform and finite difference techniques, respectively. They found that the magnetic field resists the fluid flow due to the Lorentz force. In contrast, the thermal radiation and Newtonian heating parameters lead to the enhancement of velocity and temperature fields.
The impact of the Cattaneo-Christov heat flux in Jeffrey's fluid flow with homogeneous-heterogeneous reactions was analyzed by Hayat et al. [23]. The series results they obtained show that the ratio of relaxation to retardation times and Deborah number have inverse relation for the velocity profile, and temperature distribution has decreasing behavior for the Prandtl number and thermal relaxation time. They also reported that concentration decreases for larger values of strength of the homogeneous reaction parameter while it increases for the power of the heterogeneous reaction parameter. Hayat et al. [24] considered the unsteady flow and heat transfer of Jeffrey's fluid over a stretching sheet using the homotopy analysis method. They obtained exact solutions of the momentum equation and numerical solutions of the dimensionless energy equations for the steady-state case. Their results indicated that an increase in the elastic parameter of Jeffrey's fluid (Deborah number) corresponds to the rise in the velocity and the boundary layer thickness. The homotopy analysis method was used by Hayat and Mustafa [25] to analyze the effect of thermal radiation on the unsteady mixed convection flow of Jeffrey's fluid over a porous vertical stretching surface. They reported that the heat transfer rate is increased with increasing values of the Prandtl number, the radiation parameter, and the unsteadiness parameter. Muhammad et al. [26] discussed the hydromagnetic unsteady squeezing flow of Jeffrey's fluid between two parallel plates. The nonlinear governing system of ordinary differential equations was solved analytically using the homotopy analysis method and numerically using the NDsolve. They showed that an increase in the squeezing parameter enhances the velocity profile for both suction and blowing cases. The unsteady axisymmetric flow of Jeffrey's fluid between two parallel disks was investigated by Qayyum et al. [27] using the homotopy analysis method. It was found in their investigation that the velocity profile increases when porosity and squeezing parameters are increased. Venkateswara et al. [28] studied the unsteady MHD mixed convection Jeffrey's fluid flow over an inclined permeable moving plate in the presence of thermal radiation, heat generation, thermophoresis effect, and homogenous chemical reaction, subjected to variable suction. In their study, they solved the governing equations using a regular perturbation technique. The study showed that the velocity increases with an increase in the 2 Abstract and Applied Analysis Soret number in the presence of permeability. At the same time, a reverse effect in the case of the heat absorption coefficient, magnetic parameter, radiation parameter, and chemical reaction parameter was observed.
In the past few decades, the study of the flow and the thermophysical properties of nanofluids has been given much attention due to the enormous importance of these fluids, such as in microelectronics, fuel cells, pharmaceuticals processes, cancer therapy, and electronics, among others. The concept of nanofluids was first introduced by Choi [29], where he proposed the suspension of nanoparticles in a base fluid such as water, oil, and ethylene glycol. Buongiorno [30] developed a nonhomogeneous two-component equation for nanofluids. He proposed seven slip mechanisms between nanoparticles and the base fluid. He attempted to explain the increase in the thermal conductivity of nanofluids, and his model took into account the particle Brownian motion and thermophoresis. Dhanai et al. [31] presented multiple solutions of MHD boundary layer flow and heat transfer behavior of nanofluids induced by a power-law stretching/shrinking permeable sheet with viscous dissipation with the aid of the shooting method. In their study, viscous dissipation was found to be significant, whereas the Brownian motion has a negligible effect on the rate of heat transfer. The generalized diffusion effects on the Maxwell nanofluid stagnation point flow over a stretchable sheet with slip conditions and the chemical reaction were analyzed numerically by Khan et al. [8] using the fourth-order Runge-Kutta method along with the shooting technique. Their results showed that the skin friction coefficient reduces for large values of the slip parameter, but opposing behavior is noticed for the fluid relaxation parameter. The homotopy analysis method was used by Alamri et al. [32] to analyze the effects of the second-order slip-on plane Poiseuille nanofluid under the influence of the Stefan blowing in a channel. They observed that the slowing down effects of the Stefan blowing are significantly seen for velocity and temperature profiles, whereas the opposite characteristic for the case of nanoparticle concentrations is noticed. Also, an extra sensitivity in the field of velocity was observed for a secondorder slip as compared to the first-order slip.
Nadeem et al. [33] used the homotopy analysis method to solve the nonlinear differential equations governing the oblique stagnation point flow of the Casson nanofluid towards a stretching surface with heat transfer. They found that the boundary layer is formed when the surface's stretching velocity is less than the inviscid free-stream velocity, and velocity at a point decreases with the increase in a non-Newtonian (Casson) parameter of the fluid. An optimal and numerical tactic was used by Shah et al. [34] to get the solution of radiative heat and mass transfer analysis of the micropolar nanofluid flow of the Casson fluid between two rotating parallel plates with effects of the Hall current. They found that an amassed Hall impact decreases the effective conductivity, which intends to increase the velocity field, and the temperature field enhances with larger values of the Brownian motion thermophoresis effect. Recently, Nadeem and Khan [35] obtained dual solutions of MHD oblique stagnation point flow of nanofluid over an oscillatory stretching/shrinking sheet using the bvp4c package in MATLAB. They showed that dual solutions exist in both the stretching case and the shrinking case. Also, the lower solution branch they obtained shows singular behavior for the skin friction coefficient in the shrinking domain. In contrast, the upper solution branch has smooth behavior in both the stretching and shrinking domains. Other recent studies of nanofluid flows include those by Yousif et al. [36], Sadiq et al. [37], Besthapu et al. [38], Kamal et al. [39], Nasir et al. [40], Nayak et al. [41], Salawu and Ogunseye [42], Kumar and Srinivas [43], and Awais et al. [44], among others.
As far as the authors are aware, the study of Jeffrey's nanofluid has not been given much attention. In this study, the fluid flow, heat, and mass transfer in Jeffrey's fluid containing suspensions of nanoparticles are studied with combined effects of various parameters entering the flow problem. This work is aimed at studying the effects of heat generation, thermal radiation, chemical reaction, and some other parameters discussed in the problem on natural convection viscoelastic Jeffrey's nanofluid flow from a vertical permeable flat plate. The traditional model of Hussain et al. [45] is revised into Jeffrey's nanofluid model. The flow is characterized by thermal radiation, suction/injection, chemical reaction Brownian motion, and a thermophoretic force. Much emphasis is laid on the nanoparticle concentration boundary conditions, which is made to be more realistic due to the incorporation of the Brownian and thermophoresis diffusion coefficients, and their effect on the fluid model is studied. Another aim of the study is to apply, for the first time, the multidomain bivariate spectral quasilinearization method (MD-BSQLM) to solve coupled nonlinear systems of partial differential equations (PDEs) describing Jeffrey's nanofluid-type equations. The multidomain bivariate spectral quasilinearization method linearizes the nonlinear governing system of PDEs using the Newton-Raphson-based quasilinearization method of Bellman and Kalaba [46]. The resulting equations are then integrated into multiple subintervals using the Chebyshev spectral collocation method with the Lagrange interpolation polynomials as basis functions. The Chebyshev spectral collocation method with the Lagrange interpolation polynomials are applied to the linearized nonlinear systems of partial differential equations independently in both space and time directions. These useful features of the MD-BSQLM enable the approach to yield a very accurate solution. The method has a much better region of convergence for the approximate solution when compared to other Chebyshev spectral collocationbased methods such as bivariate spectral homotopy analysis method [47], bivariate spectral quasilinearization [48], and bivariate spectral relaxation method [49], among others. This study sought, among other things, to check the accuracy and robustness of the MD-BSQLM scheme in finding solutions to this class of problems with significant complexities. Accuracy of the numerical method used is established by computing and analyzing solutions of the momentum, heat, and mass equation. Solutions found using the multidomain bivariate spectral quasilinearization method are compared with solutions existing in the 3 Abstract and Applied Analysis literature, and a good agreement is observed. The effect of different embedded physical parameters on the velocity, temperature, and concentration fields is discussed in tabular and graphical forms. The results illustrate the different behavior that occurs when these parameters are varied. The present study finds applications in polymeric manufacturing processes, heat exchanger technology nuclear waste simulations, nuclear engineering, thermal fabrication of paint sprays, and low-density polymeric materials in the process engineering industry.

Mathematical Formulation of the Problem
A steady two-dimensional natural convection boundary layer flow of an incompressible viscoelastic Jeffrey's nanofluid from a vertical permeable flat plate, as shown in Figure 1, is investigated. The fluid is maintained at the same temperature and concentration. The temperature and concentration of the fluid are raised simultaneously. The acceleration due to gravity g acts vertically downward. The surface of the fluid is held at a variable temperature and concentration proportional to the power of the distance along the surface; i.e., T w ðxÞ = T ∞ + Bd 1xn and C w ðxÞ = C∞+Bd 2 x n , where B, d 1 , d 2 are constants and n is the power-law exponent. As highlighted by Bird et al. [50], the physical characteristics of certain polymers are accurately captured by Jeffrey's model. The Cauchy stress tensor, S, of Jeffrey's viscoelastic nanofluid is given by [2] as where a dot above a quantity signifies the material time derivative, _ γ is the shear rate, μ is the dynamic viscosity, λ is the ratio of relaxation to retardation times, and λ 1 is the retarda-tion time. The shear rate and gradient of the shear rate are defined in terms of velocity vector V by [2] as Under the Boussinesq boundary layer approximations, the flow is governed by the following equations (see [2,45,51]) as In the above equations, u and v are the velocity components in x and y directions, respectively, ν = ðμ/ρÞ is the kinematic coefficient of viscosity of the fluid, g is the acceleration due to gravity, β and β * are the thermal expansion and concentration expansion coefficients, respectively, T is the fluid temperature, C is the fluid concentration, T ∞ is the ambient fluid temperature, C ∞ is the ambient fluid concentration, k 0 is the thermal diffusivity ratio, ρ is the density of the fluid, c p is the specific heat at constant pressure of the fluid, q r is the radiative heat flux, Q 0 is the heat generation constant, D m is the mass diffusivity, c s is the concentration susceptibility, τ = ðρcÞ p / ðρcÞ f is the ratio of the heat capacity of the nanoparticle material and the heat capacity of the fluid, D B is the Brownian diffusion coefficient, D T is the thermophoretic diffusion coefficient, T m is the mean fluid temperature, and K r is the chemical reaction coefficient. The radiative heat flux q r is defined using the Rosseland approximation [52] as where σ * and k * are the Stefan-Boltzmann constant and mean absorption coefficient, respectively. Assuming that the temperature differences within the flow are sufficiently Figure 1: The flow coordinate system and the flow configuration. 4 Abstract and Applied Analysis small, T 4 may be approximated in the Taylor series form about T ∞ after ignoring higher-order terms as The associated boundary conditions are where k ∞ is the thermal conductivity, h f is the convective heat transfer coefficient, and V w is the transpiration velocity of the fluid. V w > 0 denotes suction or withdrawal, i.e., mass flux removal from the boundary layer through the permeable surface wall into the permeable surface, and V w < 0 stands for injection or blowing of fluid through the permeable surface. In this present investigation, the case of suction or injection will only be considered rather than the blowing case, and therefore, V w is taken to be positive throughout this study. The following dimensionless transformations are then introduced into equations (4)-(9): In equation (10), Gr x is the local Grashof number, ξ is the transpiration parameter depending on the transpiration velocity V w , x is the axial variable, η is the pseudosimilarity variable, f , θ, φ are the dimensionless velocity, dimensionless temperature, and dimensionless concentration of the fluid, respectively, and ψ is the stream function defined as which satisfies the continuity equation (1). Substituting the transformations given in equation (10) into equations (4)-(9) gives the following set of the nonsimilarity system of partial differential equations which are expressed in a dimensionless form as subject to the boundary conditions: where the prime denotes differentiation with respect to η, x /x 2 is the Deborah number, λ is the ratio of relaxation to retardation time, NR = 16σ * T 3 ∞ /3k * vρc p k 0 is the radiation parameter, Pr = vρc p /k 0 is the Prandtl number, n is the power-law exponent, Abstract and Applied Analysis number, δ = K r x 2 /vGr 1/2 x is the chemical reaction parameter, f w = ð4/3 + nÞðxV w /νGr 1/4 x Þ is the suction or injection parameter, and Bi = xh f /k * Gr 1/4 x is the Biot number. The physical quantities of importance are the skin friction, the Nusselt number, and the Sherwood number which may be obtained from

Multidomain Bivariate Spectral Quasilinearization Method
In this section, a description of the multidomain bivariate spectral quasilinearization method (MD-BSQLM) for solving the set of coupled nonlinear partial differential equations (17)-(19) is given. The basic idea of the approach is to decompose solutions over a large interval into smaller subintervals, then solve the system of equations in the small subintervals, and take the solution at the end of each of the subintervals to get the resulting solution. The multidomain approach is applied in the ξ direction only. We first linearize equations (17)-(19) using the quasilinearization (QLM) as described in [46,53]. This technique uses the Taylor series expansion to linearize nonlinear differential equations. In the linearization technique, we make the assumption that the difference between the value of the unknown function at the current iteration level denoted by r + 1 and the value at the previous iteration level denoted by r is small. Applying the QLM on (17)- (19) gives the following: a 9,r = 1, Abstract and Applied Analysis Sc ∂f r ∂ξ , Scf r ′, Sc ∂ϕ ∂ξ , Scϕ r ′ , Sc Now, let ξ ∈ Ω, where Ω ∈ ½0, T, and the domain Ω is decomposed into p nonoverlapping intervals as The partial differential equations are solved independently at each of the p subintervals. Once the solution at the first subinterval has been computed, the new solutions at the subsequent mth interval are computed using the solution at the right-hand boundary of the m − 1th interval as an initial solution. In the mth subinterval, we solve subject to the boundary conditions: A suitable initial condition to start the multidomain iteration scheme in the first subinterval is the one that 7 Abstract and Applied Analysis satisfies the boundary conditions (20). The initial condition at the subsequent subintervals is given by the continuity conditions: The physical domains in η and ξ are first transformed to the computational domain ðx, tÞ ∈ ½−1, 1 × ½−1, 1 at each subinterval using the linear transformation: where L x is a number large enough to approximate conditions at infinity in η. The collocation points are the Chebyshev-Gauss-Lobatto nodes defined in [54,55] by where ðN x + 1Þ and ðN t + 1Þ are the total number of collocation points in η and ξ directions, respectively. Suppose that the solutions f , θ, and φ can be approximated at each subinterval by a bivariate Lagrange interpolation polynomial to the form: where the functions L p ðxÞ and L q ðtÞ are the Lagrange cardinal polynomials defined as with The first spatial derivatives of f , θ, and φ with respect to η at the Chebyshev-Gauss-Lobatto points ðx i , t j Þ for i = 0, 1, 2, ⋯, N x are evaluated as whereD = L x D/2 is the standard first derivative Chebyshev differentiation matrix of size ðN x + 1Þ × ðN x + 1Þ as 8 Abstract and Applied Analysis defined in Trefethen [54]. The vectors F and the superscript T here denotes matrix transpose. The nth-order derivatives of f , θ, and φ with respect to η are approximated using the matrix product as The spatial derivatives of f , θ, and φ are evaluated at the Chebyshev-Gauss-Lobatto points (x i , t j ) for j = 0, 1, 2 , ⋯, N t as whered j,q = ðξ m − ξ m−1 Þ/2d j,q , j, q = 0, 1, 2, N t are the entries of the standard first-order Chebyshev differentiation matrix in the mth subinterval.  Table 1: Comparison of the MD-BSQLM approximate solutions of −θ′ð0, ξÞ and −φ′ð0, ξÞ, against those of Ref. [45] for different values of n and ξ when Pr = 0:72, Sc = 0:94, N = 1/3, λ = δ = He = NR = De = Df = f w = 0, and Sr = 0 in the absence of Nt, Nb, and Bi.
Ref. [ Table 3: Table of  N Sc Pr Noting that the solution at the time level j = N t of each subinterval is given by the solution at the previous level and taking i and j as dummy indices, equations (39)- (41) can be written as In a more compact format, equations (42)-(44) can be written as Table 4: Computed values of the skin friction coefficient, local Nusselt number, and local Sherwood number for different values of Sr, Df , NR, and δ when De = 0:1, λ = 0:2, n = 0:5, Bi = 0:8, N = 0:5, ξ = 2, Pr = Sc = 1, Nb = 0:3, Nt = 0:2, and f w = 0:8.

Df
Sr The boundary conditions given in equation (28) when evaluated at the Chebyshev-Gauss-Lobatto collocation points give

Results and Discussion
In this section, we present the numerical solutions of the coupled system of nonlinear differential equations that model natural convection viscoelastic Jeffrey's nanofluid flow from a vertical permeable flat plate with thermal radiation. The modeled equations have been solved using the multidomain bivariate spectral quasilinearization method. To validate our numerical method and to show that our approach is suitable for the problem under investigation, a comparison with results in the published literature was made, and our results are in excellent agreement with these published results. We further present and examine the impact of a wide range of embedded pertinent thermophysical parameters on the flow properties, the skin friction coefficient, the Nusselt number, and the Sherwood number. Comprehensive results are obtained and are presented in Tables 1-4 and Figures 2-15. It is worth noting that the range of the parameter values used in this study is between ½0, 4. These values have been chosen based on the works of existing literature similar to this study. Table 1 shows a comparison of the Nusselt number and Sherwood number with the results of Hussain et al. [45] for different values of n and ξ. A good agreement is observed. It is seen from the table that an increase in the Nusselt number and the Sherwood number is observed as n and ξ values increase. This gives confidence that the present numerical results are accurate.
To further ascertain the accuracy of the present numerical method, we have compared the skin friction coefficient, Nusselt number, and Sherwood number for different values of N, Sc, Pr, and f w when ξ = 3 with the published result of Gaffar 13 Abstract and Applied Analysis et al. [2]. The comparison is shown in Table 2. From the table, it is seen that an excellent agreement between the present numerical results and the published results is achieved, thus validating the accuracy of the current numerical method. Also, it can be observed from the table that an increase in the values of N increases the skin friction coefficient, heat transfer rate, and mass transfer rate. The heat transfer rate increases significantly with increasing values of Pr, while a reduction in the skin friction coefficient and the mass transfer rate is observed for an increase in Pr. Also, increasing f w is seen to reduce the skin friction coefficient and heat transfer rate, while the mass transfer rate is enhanced for an increase in f w . Lastly, a decrease in the skin friction coefficient is seen with increasing values of Sc. A reduction in the heat transfer rate is observed with increasing amounts of Sc, while the mass transfer is seen to increase significantly. Table 3 shows the skin friction coefficient, Nusselt number, and Sherwood number for various values of N, Sc, Pr, f w , and ξ. It is noticed in the table that as we increase the values of N, the skin friction coefficient and the Sherwood number increase while the Nusselt number decreases. Also, as the values of Pr increase, the Nusselt number increases while the skin friction coefficient and the Sherwood number decrease. For higher values of Pr, velocity decreases, thereby reducing the skin friction coefficient. The reduction is due to the increase in the viscosity, which in turn enhances the momentum boundary layer thickness. Increasing the values of f w is seen to reduce the skin friction coefficient and mass transfer rate, while the heat transfer rate is increased. Furthermore, a significant decrease in the skin friction coefficient is observed with increasing Sc. A slight decrease in the Nusselt number is observed with increasing values of Sc 14 Abstract and Applied Analysis while the Sherwood number is seen to increase. Increasing Sc results in a decrease in species mass diffusivity. Table 4 displays the computed values of the skin friction coefficient, the heat transfer coefficient, and the mass transfer coefficient for different values of Df , Sr, NR, δ, and He. We note that an increase in Df reduces both the skin friction coefficient and the Sherwood number, whereas the Nusselt number increases when Df is increased. The skin friction coefficient and the Sherwood number are enhanced by increasing the values of NR and Sr, while the Nusselt number decreases when increasing the values of NR and Sr. The Nusselt number is slightly increased by increasing δ, but the opposite trend is observed in the case of the skin friction coefficient and the Sherwood number. Again, the skin friction coefficient and the Sherwood number increase, whereas the Nusselt number decreases when increasing the values of He.  (17); hence, the parameter De exerts a strong impact on the shearing characteristic of the flow. This will, in turn, correspond to a decrease in the momentum boundary layer thickness. A similar trend was observed in [2]. 15 Abstract and Applied Analysis relaxation to retardation time). Here, the velocity profile is significantly increased with increasing values of λ, especially that close to the surface of the wall, whereas temperature and concentration profiles slightly reduce with increasing λ. The parameter λ is expected to expend a reasonable impact on the flow properties as increasing relaxation times (or decreasing retardation times) helps in the development of momentum in the boundary layer but reduces both thermal and mass diffusion. This flow pattern with a similar parameter effect is observed in [2]. Figures 4(a)-4(c) illustrate the variation of the velocity, temperature, and concentration profiles for different values of n. We note that there is a considerable decrease in the velocity profile with n. The same trend is obtainable for the temperature and concentration profiles as they are seen to reduce with an increase in n, as presented in Figures 4(b) and 4(c). For n > 0, the wall temperature upsurges with distance from the leading edge; for n < 0, the wall temperature reduces; and for n = 0, the wall temperature is seen to be isothermal. Increasing the values of n results in a reduction in boundary

16
Abstract and Applied Analysis layer flow and a corresponding increase in momentum boundary layer thickness and a deceleration in thermal boundary layer thickness. This trend was observed in [2]. Figures 5(a)-5(c) present the influence of N on the velocity, temperature, and concentration profiles. An increase in N is observed to upsurge the velocity profile, while a decrease in both temperature and concentration profiles is seen for increased values of N. This same trend was observed in the fluid model studied by [2].
Figures 6(a)-6(c) give a visual display of the velocity, temperature, and concentration profiles for different values of f w . It can be seen from Figures 6(a)-6(c) that there is a decrease in the velocity, temperature, and concentration pro-files with an increase in f w values. The parameter f w (suction) causes the boundary layer to be closely attached to the wall, thereby reducing the momentum boundary layer thickness. The thickening of the momentum boundary layer simultaneously obstructs heat diffusion, which consequently leads to a decrease in the temperature profile. This trend was observed in [2,56].   Figures 7(b) and 7(c) that the temperature and concentration profiles are reduced by increasing ξ. In the study by [2], similar profile patterns were obtained and explained extensively. Figures 8(a) and 8(b) show the effect of Df on the temperature and concentration profiles. Increasing the Df parameter leads to a decrease in the temperature and concentration profiles. The Dufour effect is the heat transfer induced by volume fraction gradients and significant due to the density dif-ference in the flow regime. The thermal energy flux is seen to be affected by the concentration gradients because of the reduction in the temperature profile. The temperature gradient effect is reduced, thereby leading to the cooling of the boundary layer region. In effect, this means that the thermal and mass boundary layer thickness is reduced as the parameter is enhanced.

Abstract and Applied Analysis
Bi raises both the temperature and concentration profiles. An increase in Bi yields more vigorous convection, which in turn causes higher surface temperature and an increase in wall temperature, thereby leading to thicker related layers. The effect of Nt on the temperature and concentration profiles is shown in Figures 10(a) and 10(b). As Nt increases, the temperature profile decelerates, whereas the concentration profile upsurges. The thermophoresis parameter is defined as the ratio of thermophoretic diffusion due to the temperature gradient to the momentum diffusion in the nanofluid. So, a higher thermophoresis parameter implies a stronger thermophoretic force. As a result, thermophoretic force drags enormous nanoparticles through greater diffusion from the hot surface towards the ambient surface, thereby increasing the   fluid temperature and diminishing the related layer thickness. Also, higher thermophoresis values suggest a higher concentration of the nanoparticles in the fluid, hence increasing the concentration boundary layer thickness. Figures 11(a) and 11(b) show the effect of Nb on the temperature and the concentration profiles. It was observed that the temperature profile increases with Nb while the concentration profile decreases with Nb. The concentration profile is seen to decrease with increasing Nb near the wall up to a certain value of η, while beyond this point, a reverse trend is observed. Figure 12 depicts the temperature profile for various He values. An increase in He corresponds to an enhancement in temperature and boundary layer thickness. Physically, the parameter He releases heat energy to the fluid flow; this energy released enhances fluid temperature and thus leads to an increase in the thickness of the thermal boundary layer. In the investigation carried out by [51], a similar temperature profile pattern was reported. Figure 13 clearly shows that both the thermal boundary layer thickness and the temperature profile increase owing to the rise in the values of NR. A similar trend was reported in [57].

Abstract and Applied Analysis
The effect of Sr on the concentration profile is observed in Figure 14. The parameter Sr is the ratio of the thermal diffusion coefficient and diffusion coefficient. This ratio acts as the mass flux produced by temperature gradients. Therefore, as Sr increases, the concentration profile significantly increases. This implies that the temperature gradients render a mass flux, which increases the concentration profile, thus enhancing the nanoparticle concentration boundary layer thickness. Figure 15 shows the impact of δ on the concentration profiles. It is observed that increasing δ reduces the concentration profile. This suggests that higher values of δ lead to a decrease in the chemical molecular diffusivity, thereby causing a reduction in the concentration of the diffusing species, and a decrease in mass diffusion, hence leading to the thinning of the concentration boundary layer.

Conclusion
In this paper, we have investigated the combined effects of the Soret and Dufour numbers, thermal radiation, heat generation, and chemical reaction on natural convection viscoelastic Jeffrey's nanofluid flow from a vertical permeable flat plate. The governing nonlinear differential equations were solved using the multidomain bivariate spectral quasilinearization method. A comprehensive assessment of different thermophysical quantities is discussed graphically and in a tabular form. The method employed in this study has demonstrated its ability to solve nonlinear boundary layer equations efficiently and shows an excellent promise in simulating transport phenomena in other non-Newtonian fluids. In terms of the physical problem, the study has shown, among other things, that by increasing the Deborah number, power-law exponent, transpiration parameter, and suction parameter, we reduce the fluid velocity profile. Both the fluid temperature and thermal boundary layer thickness are found to increase with the heat generation parameter, Brownian motion parameter, Biot number, and thermal radiation parameter while they decreased with an increase in the thermophoresis parameter and Dufour number. The concentration profile reduces with the chemical reaction and Brownian motion parameters, whereas the profile is enhanced with the thermophoresis parameter and Soret and Biot numbers. These observations are parallel with earlier studies in the literature.

Data Availability
No data were used to support this study.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.