Numerical Solution of MHD Viscoelastic Nanofluid Flow over a Stretching Sheet with Partial Slip and Heat Source/Sink

Abstract

We analyze the effect of velocity slip boundary condition on the flow and heat transfer of non-Newtonian nanofluid over a stretching sheet with a heat source/sink, under the action of a uniform magnetic field, orientated normally to the plate. The Brownian motion and thermophoresis effects are also considered. The boundary layer equations governed by the partial differential equations are transformed into a set of ordinary differential equations with the help of local similarity transformations. The differential equations are solved by the variational finite element method (FEM). We have examined the effects of different controlling parameters, namely, the Brownian motion parameter, the thermophoresis parameter, uniform magnetic field, viscoelastic parameter, Prandtl number, heat source/sink parameter, Lewis number, and the slip parameter on the flow field and heat transfer characteristics. Graphical display of the numerical examination is performed to illustrate the influence of various flow parameters on the velocity, temperature, concentration, and Nusselt and Sherwood numbers distributions. The present study has many applications in coating and suspensions, cooling of metallic plate, paper production, heat exchangers technology, and materials processing exploiting.

1. Introduction

Due to enormous industrial, transportation, electronics, biomedical applications, such as in advanced nuclear systems, cylindrical heat pipes, automobiles, fuel cells, drug delivery, biological sensors, and hybrid-powered engines, the convective heat transfer in nanofluids has become a topic of great interest. Nanofluids are engineered by suspending nanoparticles with average sizes of 1–100 nm. Choi et al. [1] and Masuda et al. [2] have shown that a very small amount of nanoparticles (usually less than 5%), when dispersed uniformly and suspended stably in base fluids, can provide dramatic improvements in the thermal conductivity and in the heat transfer coefficient of the base fluid. The term nanofluids (nanoparticle fluid suspensions) was coined by Choi [3] in 1995 to describe this new class of nanotechnology-based heat transfer fluids that exhibit thermal properties superior to those of their base fluids or conventional particle fluid suspensions. The nanoparticles are typically made of oxides such as alumina, silica, titania and copper oxide, carbides, and metals such as copper and gold. Carbon nanotubes and diamond nanoparticles have also been used in nanofluids. The base fluid is usually a conventional heat transfer fluid, such as oil, water, and ethylene glycol. Other base fluids are biofluids, polymer solutions, and some lubricants.

A comprehensive survey of convective transport in nanofluids was made by Buongiorno [4]. He developed a nonhomogeneous equilibrium model for convective transport to explain the enhanced heat transfer characteristics of nanofluids, and this abnormal increase in the thermal conductivity occurs due to the presence of two main velocity-slip effects, namely, the Brownian diffusion and the thermophoretic diffusion of the nanoparticles. In a recent paper, Khan and Pop [5] have used the model of Kuznetsov and Nield [6] to study the fundamental work on the boundary layer flow of nanofluid over a stretching sheet. Makinde and Aziz [7] extended the work of Khan and Pop [5] for convective boundary conditions.

It is now a well-accepted fact that many fluids of industrial and geophysical importance are non-Newtonian. Due to much attention in many industrial applications, such as the extrusion of plastic sheets, fabrication of adhesive tapes, glass-fiber production, metal spinning, and drawing of paper films, the research on boundary layer behaviour of a viscoelastic fluid over a continuously stretching surface keeps going where the velocity of a stretching surface is assumed to be linearly proportional to the distance from a fixed origin. McCormack and Crane [8] have provided comprehensive discussion on boundary layer flow caused by stretching of an elastic flat sheet moving in its own plane with a velocity varying linearly with distance. Several researchers, namely, P. S. Gupta and A. S. Gupta [9], Dutta et al. [10], and Chen and Char [11] extended the work of McCormack and Crane [8] by including the effects of heat and mass transfer under different situations. Later on, Rajagopal et al. [12] and Chang [13] presented an analysis on flow of viscoelastic fluid over a stretching sheet. The above sources all utilize the no-slip condition. On the other hand, in certain circumstances, the partial slip between the fluid and the moving surface may occur in situations when the fluid is particulate such as emulsions, suspensions, foam and polymer solutions. In these cases, the proper boundary condition is replaced by Navier’s condition, where the amount of relative slip is proportional to local shear stress. Wang [14] discussed the partial slip effects on the planar stretching flow. Of late, Noghrehabadi et al. [15] investigated the development of the slip effects on the boundary layer flow and heat transfer over a stretching sheet.

In some practical problems such as the magnetohydrodynamic (MHD) generators, pumps, enhanced oil recovery, thermal insulators, electronic packages, and cooling of nuclear reactor, the flow of electrically conducting fluid occurs in the presence of a transverse magnetic field. Ece [16] analyzed the natural convection boundary layer flow in the presence of a transverse magnetic field with mixed thermal boundary conditions over a vertical down-pointing cone. The MHD boundary layer flow of an incompressible and electrically conducting viscoelastic fluid past a linear stretching sheet was studied by Subhas Abel et al. [17]. Kumari and Nath [18] have studied the unsteady magnetohydrodynamic flow and heat transfer of a Newtonian fluid caused by an impulsively stretched plane surface in two lateral directions. The momentum and heat transfer characteristics of the boundary layers of an incompressible electrically conducting fluid flow of a viscoelastic fluid over a stretching sheet are investigated by Prasad et al. [19]. Recently, Hamad et al. [20] reported a similarity solution for MHD free convection heat generation flow over a vertical semi-infinite flat plate in the case of nanofluids.

A study of utilizing heat source or sink in moving fluids assumes a greater significance in all situations which deal with exothermic or endothermic chemical reaction and those concerned with dissociating fluids. For physical situations, the average behaviour of heat generation or absorption can be expressed by some simple mathematical models because its exact modelling is quite difficult. Heat generation or absorption has been assumed to be constant, space dependent, or temperature-dependent. Sparrow and Cess [23] investigated the steady stagnation point flow and heat transfer in the presence of temperature dependent heat absorption. Later, Azim et al. [24] discussed the effect of viscous Joule heating on MHD-conjugate heat transfer for a vertical flat plate in the presence of heat generation. One of the latest works is the study of the heat transfer characteristic in the mixed convection flow of a nanofluid along a vertical plate with heat source/sink studied by Rana and Bhargava [25].

In real situations in nanofluids, the base fluid does not satisfy the properties of Newtonian fluids; hence, it is more justified to consider them as viscoelastic fluids; for example, ethylene glycol-Al₂O₃, ethylene glycol-CuO, and ethylene glycol-ZnO are some examples of viscoelastic nanofluids. In the present paper, the base fluid is taken as second-grade fluid. To the best of our knowledge, no studies have far been investigated to analyze the partial slip effect on the MHD boundary layer flow of viscoelastic nanofluid over a stretching sheet in the presence of heat source/sink. The objective of the present paper is therefore to extend the work of Noghrehabadi et al. [15] by taking base fluid as second-grade fluid. The finite element method (FEM) is used in obtaining the numerical solution of the obtained equations, which is describing the problem, after similarity transformation. A similarity solution is presented and used to predict the heat and mass transfer characteristics of the flow. The effects of the embedded flow controlling parameters on the fluid velocity, temperature, nanoparticle concentration, heat transfer rate, and the nanoparticle volume fraction rate have been demonstrated graphically and discussed. A comparative study is also presented.

2. Formulation

Consider two-dimensional, steady, incompressible, and laminar flow of non-Newtonian nanofluid past a stretching sheet in a quiescent fluid. The velocity of the stretching sheet is 𝑢_𝑤=𝑈=𝑐𝑥, where 𝑐>0 is the constant acceleration parameter. In order to get the effect of temperature difference between the surface and the ambient fluid, we consider temperature-dependent heat source/sink in the flow. The x-axis is taken along the plate in the vertically upward direction and the y-axis is taken normally to the plate. A transverse magnetic field of strength 𝐵_𝑜 is applied parallel to the y-axis as shown in Figure 1. The surface of plate is maintained at uniform temperature 𝑇_𝑤 and concentration 𝐶_𝑤, and these values are assumed to be greater than the ambient temperature 𝑇_∞ and concentration 𝐶_∞. Moreover, it is assumed that both the fluid phase and nanoparticles are in thermal equilibrium state. The thermophysical properties of the nanofluid are assumed to be constant. The pressure gradient and external forces are neglected. Consider 𝜕𝑢+𝜕𝑥𝜕𝑣𝑢𝜕𝑦=0,𝜕𝑢𝜕𝑥+𝑣𝜕𝑢𝜕𝜕𝑦=𝜐²𝑢𝜕𝑦²+𝛼₁𝜌_𝑓𝜕𝑢𝜕𝜕𝑥²𝑢𝜕𝑦²𝜕+𝑢³𝑢𝜕𝑥𝜕𝑦²+𝜕𝑢𝜕𝜕𝑦²𝑣𝜕𝑦²𝜕+𝑣³𝑢𝜕𝑦³−𝜎𝜌_𝑓𝐵²_𝑜𝑢𝑢,𝜕𝑇𝜕𝑥+𝑣𝜕𝑇𝜕𝑦=𝛼_𝑚∇²𝑇+𝜏𝐷_𝐵𝜕𝐶⋅𝜕𝑦𝜕𝑇+𝐷𝜕𝑦𝑇^𝑇_∞𝜕𝑇𝜕𝑦²+𝑄_𝑜𝜌𝐶_𝑓𝑇−𝑇_∞,𝑢𝜕𝐶𝜕𝑥+𝑣𝜕𝐶𝜕𝑦=𝐷_𝐵𝜕²𝐶𝜕𝑦²+𝐷_𝑇𝑇_∞𝜕²𝑇𝜕𝑦², (1) where 𝑢 and 𝑣 are the velocity components along the 𝑥 and 𝑦 directions, respectively, 𝑝 is the pressure, 𝜌_𝑓 is the density of base fluid, 𝜌_𝑝 is the nanoparticle density, 𝜇 is the absolute viscosity of the base fluid, 𝜐 is the kinematic viscosity of the base fluid, 𝜎 is the electrical conductivity of the base fluid, 𝛼₁ is the material fluid parameter, 𝐵_𝑜 is the strength of magnetic field, 𝑇 is the fluid temperature, 𝛼_𝑚 is the thermal diffusivity, 𝜏(=(𝜌𝐶)_𝑝/(𝜌𝐶)_𝑓) is the ratio of effective heat capacity of the nanoparticle material to heat capacity of the fluid, 𝐶 is the nanoparticle volume fraction, 𝐷_𝐵 and 𝐷_𝑇 are the Brownian diffusion coefficient and the thermophoresis diffusion coefficient, 𝑇_∞ is the free stream temperature, 𝐶_𝑝 is the specific heat at constant pressure, and 𝑔, 𝑘 are the acceleration due to gravity and the thermal conductivity of the fluid, respectively. The term 𝑄_𝑜(𝑇−𝑇_∞) is assumed to be the amount of heat generated or absorbed per unit volume 𝑄_𝑜 as a coefficient constant, which may take on either positive or negative value. When the wall temperature 𝑇 exceeds the free stream temperature 𝑇_∞, the source term 𝑄_𝑜>0 and heat sink when 𝑄_𝑜<0.

Figure 1 

Physical model and coordinate system.

The boundary conditions for the velocity, temperature, and concentration fields are given as follows:𝑢=𝑈+𝜅𝜐𝜕𝑢𝜕𝑦,𝑣=0,𝑇=𝑇_𝑤,𝐶=𝐶_𝑤,at𝑦=0,𝑢=0,𝑇=𝑇_∞,𝐶=𝐶_∞as𝑦⟶∞. (2a)(2b)To transform the governing equations into a set of similarity equations, the following dimensionless parameters are introduced: 𝜂=𝑐𝜐𝑦,𝑢=𝑐𝑥𝑓^′𝜂,𝑣=−√𝑐𝜐𝑓𝜂,𝜃𝜂=𝑇−𝑇_∞𝑇_𝑤−𝑇_∞,𝜙𝜂=𝐶−𝐶_∞𝐶_𝑤−𝐶_∞. (3) The transformed momentum, energy, and concentration equations together with the boundary conditions given by (1), (2a), and (2b) can be written as 𝑓^′′′+𝑓𝑓^′′−𝑓^′2−𝛼𝑓^′′2−2𝑓^′𝑓^′′′+𝑓𝑓^𝑖𝑣−𝑀𝑓^′1=0,𝜃Pr^′′+𝑓𝜃^′+Nb𝜃^′𝜙^′+Nt𝜃^′2𝜙+𝑄𝜃=0,_′′+Le𝑓𝜙^′+Nt𝜃Nb^′′=0. (4) The transformed boundary conditions are𝑓(0)=0,𝑓^′(0)=1+𝐾𝑓^′′(0),𝜃(0)(0)𝑓=1,𝜙=1at𝜂=0,_′(∞)=0,𝜃(∞)=0,𝜙(∞)=0,as𝜂⟶∞, (5a)(5b)where primes denote differentiation with respect to 𝜂 and the seven parameters appearing in (4) are defined as follows: 𝜐Pr=𝛼_𝑚𝜐,Le=𝐷_𝐵𝜎,𝑀=𝜌_𝑓𝐵²_𝑜𝑏𝛼,𝛼=¹𝑏𝜇,𝑄𝑄=^𝑜𝑏𝜌𝐶_𝑓,Nb=𝜌𝐶_𝑝𝐷_𝐵𝐶_𝑤−𝐶_∞𝜌𝐶_𝑓𝜐,Nt=𝜌𝐶_𝑝𝐷_𝑇𝑇_𝑤−𝑇_∞𝜌𝐶_𝑓𝑇_∞𝜐. (6) In (6), Pr,Le,𝑀,𝛼,𝑄,Nb, and Nt denote the Prandtl number, the Lewis number, the magnetic field strength parameter, the viscoelastic parameter, the heat source/sink parameter, the Brownian motion parameter, and the thermophoresis parameter, respectively.

The physical quantities of interest are the local heat flux Nu and the local mass diffusion flux Sh, which are defined as Nu=𝑥𝑞_𝑤𝑘𝑇_𝑤−𝑇_∞,Sh=𝑥𝑞_𝑚𝐷_𝐵𝐶_𝑤−𝐶_∞, (7) where 𝑞_𝑤 and 𝑞_𝑚 are the heat flux and the mass flux at the surface, respectively, which are given by 𝑞_𝑤=−𝑘𝜕𝑇𝜕𝑦_𝑦=0,𝑞_𝑚=−𝐷_𝐵𝜕𝐶𝜕𝑦_𝑦=0. (8) Using (3) in (7), one can obtain Re_𝑥^−1/2Nu_𝑥=−𝜃^′(0),Re_𝑥^−1/2Sh_𝑥=−𝜙^′(0), (9) where Re_𝑥=𝑢_𝑤(𝑥)𝑥/𝜐 is the local Reynolds number based on the stretching velocity 𝑢_𝑤(𝑥). Kuznetsov and Nield [6] referred to Re_𝑥^−1/2Nu_𝑥 and Re_𝑥^−1/2Sh_𝑥 as the reduced Nusselt number Nur=−𝜃^′(0) and reduced Sherwood number Shr=−𝜙^′(0), respectively. The analytical solutions of the set of partial differential equations given by (4) are generally intractable because these equations are highly nonlinear. The variational finite element method (FEM) is used to obtain numerical solutions. These techniques have been used very successfully in nonlinear magnetofluid dynamics.

3. Method of Solution

3.1. Finite Element Method

The finite element method (FEM) is a numerical and computer-based technique of solving a variety of practical engineering problems that arise in different fields such as, in heat transfer, fluid mechanics [26], chemical processing [27], rigid body dynamics [28], solid mechanics [29], and many other fields. It is recognized by developers and users as one of the most powerful numerical analysis tools ever devised to analyze complex problems of engineering. The sophistication of the method, its accuracy, simplicity, and computability all make it a widely used tool in the engineering modelling and design process. It has been applied to a number of physical problems, where the governing differential equations are solved by transforming them into a matrix equation. The primary feature of FEM is its ability to describe the geometry or the media of the problem being analyzed with great flexibility. This is because the discretization of the domain of the problem is performed using highly flexible uniform or nonuniform patches or elements that can easily describe complex shapes. The method essentially consists in assuming the piecewise continuous function for the solution and obtaining the parameters of the functions in a manner that reduces the error in the solution. The steps involved in the finite element analysis are as follows.

3.1.1. Discretization of the Domain

The basic concept of the FEM is to divide the domain or region of the problem into small connected patches, called finite elements. The collection of elements is called the finite element mesh. These finite elements are connected in a nonoverlapping manner, such that they completely cover the entire space of the problem.

3.1.2. Generation of the Element Equations

(i)A typical element is isolated from the mesh and the variational formulation of the given problem is constructed over the typical element.(ii)Over an element, an approximate solution of the variational problem is supposed, and by substituting this in the system, the element equations are generated.(iii)The element matrix, which is also known as stiffness matrix, is constructed by using the element interpolation functions.

3.1.3. Assembly of the Element Equations

The algebraic equations so obtained are assembled by imposing the interelement continuity conditions. This yields a large number of algebraic equations known as the global finite element model, which governs the whole domain.

3.1.4. Imposition of the Boundary Conditions

On the assembled equations, the Dirichlet and Neumann boundary conditions are imposed.

3.1.5. Solution of Assembled Equations

The assembled equations so obtained can be solved by any of the numerical techniques, namely, Gauss elimination method, LU decomposition method, and so forth.

For the solution of system of simultaneous partial differential equations as given in (4), with the boundary conditions (5a), (5b), firstly we assume that 𝑓^′=ℎ. (10) After substituting (10), the system of (4) reduces to ℎ^′′+𝑓ℎ^′−ℎ²−𝛼ℎ^′2−2ℎℎ^′′+𝑓ℎ^′′′1−𝑀ℎ=0,𝜃Pr^′′+𝑓𝜃^′+Nb𝜃^′𝜙^′+Nt𝜃^′2𝜙+𝑄𝜃=0,_′′+Le𝑓𝜙^′+Nt𝜃Nb^′′=0, (11) and the corresponding boundary conditions are: 𝑓(0)=0,ℎ(0)=1+𝐾ℎ^′(0),𝜃(0)(0)ℎ(∞)(∞)𝜙(∞)=1,𝜙=1,at𝜂=0,⟶0,𝜃⟶0,⟶0,as𝜂⟶∞. (12a)(12b)

3.2. Variational Formulation

The variational form associated with (10) and (11) over a typical quadratic element (𝜂_𝑒,𝜂_𝑒+1) is given by ^𝜂𝑒+1_𝜂𝑒𝑊₁𝑓^′−ℎ𝑑𝜂=0,_{𝜂𝑒+1𝜂𝑒}𝑊₂ℎ^′′+𝑓ℎ^′−ℎ²−𝛼ℎ^′2−2ℎℎ^′′+𝑓ℎ^′′′−𝑀ℎ𝑑𝜂=0,_{𝜂𝑒+1𝜂𝑒}𝑊₃1𝜃Pr^′′+𝑓𝜃^′+Nb𝜃^′𝜙^′+Nt𝜃^′2+𝑄𝜃𝑑𝜂=0_{𝜂𝑒+1𝜂𝑒}𝑊₄𝜙^′′+Le𝑓𝜙^′+Nt𝜃Nb^′′𝑑𝜂=0, (13) where 𝑊₁,𝑊₂,𝑊₃, and 𝑊₄ are arbitrary test functions and may be viewed as the variations in 𝑓,ℎ,𝜃, and 𝜙, respectively.

3.3. Finite Element Formulation

The finite element model may be obtained from above equations by substituting finite element approximations of the form 𝑓=³_𝑗=1𝑓_𝑗𝜓_𝑗,ℎ=³_𝑗=1ℎ_𝑗𝜓_𝑗,𝜃=³_𝑗=1𝜃_𝑗𝜓_𝑗,𝜙=³_𝑗=1𝜙_𝑗𝜓_𝑗, (14) with 𝑊₁=𝑊₂=𝑊₃=𝑊₄=𝜓_𝑗,𝑗=1,2,3. (15) In our computations, the shape functions for a typical element (𝜂_𝑒,𝜂_𝑒+1) are taken as follows. Quadratic element: 𝜓^𝑒₁=𝜂_𝑒+1−𝜂_𝑒−2𝜂𝜂_𝑒+1−𝜂𝜂_𝑒+1−𝜂_𝑒²,𝜓^𝑒₂=4𝜂−𝜂_𝑒𝜂_𝑒+1−𝜂𝜂_𝑒+1−𝜂_𝑒²,𝜓^𝑒₃𝜂=−^𝑒+1−𝜂_𝑒−2𝜂𝜂−𝜂_𝑒𝜂_𝑒+1−𝜂_𝑒²,𝜂_𝑒≤𝜂≤𝜂_𝑒. (16) The finite element model of the equations thus formed is given by 𝐾¹¹𝐾¹²𝐾¹³𝐾¹⁴𝐾²¹𝐾²²𝐾²³𝐾²⁴𝐾³¹𝐾³²𝐾³³𝐾³⁴𝐾⁴¹𝐾⁴²𝐾⁴³𝐾⁴⁴𝑓{ℎ}{𝜃}𝜙=𝑏¹𝑏²𝑏³𝑏⁴, (17) where [𝐾^𝑚𝑛] and [𝑏^𝑚](𝑚,𝑛=1,2,3,4) are defined as 𝐾¹¹_𝑖𝑗=^𝜂𝑒+1_𝜂𝑒𝜓_𝑖𝜕𝜓_𝑗𝜕𝜂𝑑𝜂,𝐾¹²_𝑖𝑗=−^𝜂𝑒+1_𝜂𝑒𝜓_𝑖𝜓_𝑗𝐾𝑑𝜂,_13𝑖𝑗=𝐾¹⁴_𝑖𝑗=0,𝐾²¹_𝑖𝑗=𝐾²³_𝑖𝑗=𝐾²⁴_𝑖𝑗𝐾=0,_22𝑖𝑗=−^𝜂𝑒+1_𝜂𝑒𝜕𝜓_𝑖𝜕𝜂𝜕𝜓_𝑗𝜕𝜂𝑑𝜂+^𝜂𝑒+1_𝜂𝑒𝜓_𝑖𝑓𝜕𝜓_𝑗𝜕𝜂𝑑𝜂−^𝜂𝑒+1_𝜂𝑒𝜓_𝑖ℎ𝜓_𝑗𝑑𝜂−𝑀^𝜂𝑒+1_𝜂𝑒𝜓_𝑖𝜓_𝑗𝑑𝜂+𝛼2^𝜂𝑒+1_𝜂𝑒𝜓_𝑖ℎ𝜕²𝜓_𝑗𝜕𝜂²𝑑𝜂−^𝜂𝑒+1_𝜂𝑒𝜓_𝑖ℎ^′𝜕𝜓_𝑗+𝜕𝜂𝑑𝜂_{𝜂𝑒+1𝜂𝑒}𝜕𝜓_𝑖𝜕𝜂𝑓𝜕²𝜓_𝑗𝜕𝜂²𝑑𝜂+^𝜂𝑒+1_𝜂𝑒𝜓_𝑖𝑓^′𝜕²𝜓_𝑗𝜕𝜂²,𝐾𝑑𝜂_33𝑖𝑗1=−Pr^𝜂𝑒+1_𝜂𝑒𝜕𝜓_𝑖𝜕𝜂𝜕𝜓_𝑗𝜕𝜂𝑑𝜂+^𝜂𝑒+1_𝜂𝑒𝜓_𝑖𝑓𝜕𝜓_𝑗𝜕𝜂𝑑𝜂+Nb^𝜂𝑒+1_𝜂𝑒𝜓_𝑖𝜙^′𝜕𝜓_𝑗𝜕𝜂𝑑𝜂+Nt^𝜂𝑒+1_𝜂𝑒𝜓_𝑖𝜃^′𝜕𝜓_𝑗𝜕𝜂𝑑𝜂+𝑄^𝜂𝑒+1_𝜂𝑒𝜓_𝑖𝜓_𝑗𝐾𝑑𝜂,_31𝑖𝑗=𝐾³²_𝑖𝑗=𝐾³⁴_𝑖𝑗=0,𝐾⁴¹_𝑖𝑗=𝐾⁴²_𝑖𝑗𝐾=0,_43𝑖𝑗=−NtNb^𝜂𝑒+1_𝜂𝑒𝜕𝜓_𝑖𝜕𝜂𝜕𝜓_𝑗𝐾𝜕𝜂𝑑𝜂,_44𝑖𝑗=−^𝜂𝑒+1_𝜂𝑒𝜕𝜓_𝑖𝜕𝜂𝜕𝜓_𝑗𝜕𝜂𝑑𝜂+Le^𝜂𝑒+1_𝜂𝑒𝜓_𝑖𝑓𝜕𝜓_𝑗𝜕𝜂𝑑𝜂, (18) where 𝑓=³_𝑖−1𝑓_𝑖𝜓_𝑖,𝑓^′=³_𝑖=1𝑓_𝑖𝜕𝜓_𝑖,𝜕𝜂ℎ=³_𝑖−1ℎ_𝑖𝜓_𝑖,ℎ^′=³_𝑖=1ℎ_𝑖𝜕𝜓_𝑖,𝜕𝜂𝜃^′=²_𝑖=1𝜃_𝑖𝜕𝜓_𝑖,𝜕𝜂𝜙^′=²_𝑖=1𝜙_𝑖𝜕𝜓_𝑖.𝜕𝜂 (19) In one-dimensional space, linear element, quadratic element, or element of higher order can be taken. The entire flow domain is divided into 10000 quadratic elements of equal size. Each element is three-noded, and therefore the whole domain contains 20001 nodes. At each node, four functions are to be evaluated; hence, after assembly of the element equations, we obtain a system of 80004 equations which are nonlinear. Therefore, an iterative scheme must be utilized in the solution. After imposing the boundary conditions, a system of equations has been obtained which is solved by the Gauss elimination method while maintaining an accuracy of 0.00001. A convergence criterion based on the relative difference between the current and previous iterations is employed. When these differences satisfy the desired accuracy, the solution is assumed to have been converged and iterative process is terminated. The Gaussian quadrature is implemented for solving the integrations. The code of the algorithm has been executed in MATLAB running on a PC. Excellent convergence was achieved for all the results.

4. Results and Discussion

The nonlinear ordinary differential equations (4) together with the boundary conditions (5a) and (5b) are solved numerically using FEM. The numerical computations have been carried out for different values of the parameters involved, namely, viscoelastic parameter 𝛼, Prandtl number Pr, magnetic field parameter 𝑀, slip parameter 𝐾, Lewis number Le, heat source/sink parameter 𝑄, the Brownian motion parameter Nb, and thermophoresis parameter Nt. The effects of the flow controlling parameters on the dimensionless axial velocity, temperature, the nanoparticle volume fraction, and the rate of heat and mass transfer are investigated and presented graphically in Figures 2–15. The aim of the present study is to examine the variations of different quantities of parameters in which 0≤𝐾≤5,0≤𝛼≤10,0≤Pr≤70,0≤𝑀≤5,0.1≤Nt≤0.5,0.1≤Nb≤0.5,−0.5≤𝑄≤0.5, and 5≤Le≤20. The computational work is carried out by taking size of the element ∇𝜂=0.0001. The accuracy of current solutions is independent of the size of the element, thus, the results are convergence solutions.

Figure 2 

Velocity profile for various values of 𝐾 when Pr=Le=10,Nb=Nt=0.1, 𝑀=1.0,𝑄=0.05,𝛼=0.5.

Figure 3 

Temperature profile for various values of 𝐾 when Pr=Le=10,Nb=Nt=0.1,𝑀=1.0,𝑄=0.05,𝛼=0.5.

Figure 4 

Nanoparticle concentration profile for various values of 𝐾 when Pr=Le=10,Nb=Nt=0.1,𝑀=1.0,𝑄=0.05,𝛼=0.5.

Figure 5 

Velocity profile for various values of 𝑀 when Pr=Le=10,Nb=Nt=0.1,𝐾=1.0,𝑄=0.1,𝛼=0.5.

Figure 6 

Temperature profile for various values of 𝑀 when Pr=Le=10,Nb=Nt=0.1,𝐾=1.0,𝑄=0.1,𝛼=0.5.

Figure 7 

Nanoparticle concentration profile for various values of 𝑀 when Pr=Le=10,Nb=Nt=0.1,𝐾=1.0,𝑄=0.05,𝛼=0.5.

Figure 8 

Velocity profile for various values of 𝛼 when Pr=5,Le=10,Nb=Nt=0.1,𝐾=1.0,𝑄=0.1,𝑀=1.0.

Figure 9 

Temperature profile for various values of 𝛼 when Pr=5,Le=10,Nb=Nt=0.1,𝐾=1.0,𝑄=0.1,𝑀=1.0.

Figure 10 

Nanoparticle concentration profile for various values of 𝛼 when Pr=5,Le=10,Nb=Nt=0.1,𝐾=1.0,𝑄=0.1,𝑀=1.0.

Figure 11 

Temperature profile for various values of 𝑄 when Pr=Le=10,𝑀=1.0,Nb=Nt=0.1,𝐾=1.0,𝛼=0.5.

Figure 12 

Heat transfer rate against Nt for several values of 𝐾 when Pr=Le=10,Nb=0.1,𝛼=0.5,𝑀=1.0,𝑄=0.1.

Figure 13 

Mass transfer rate against Nt for several values of 𝐾 when Pr=Le=10,Nb=0.1,𝛼=0.5,𝑀=1.0,𝑄=0.1.

Figure 14 

Heat transfer rate against Nt for several values of 𝛼 when Pr=Le=10,Nb=0.1,𝑀=0.0,𝑄=0=𝐾.

Figure 15 

Mass transfer rate against Nt for several values of 𝛼 when Pr=Le=10,Nb=0.1,𝑀=0.0,𝑄=0=𝐾.

Figures 2, 3, and 4 illustrate the velocity, temperature, and concentration profiles for different values of the slip parameter 𝐾. Figure 2 demonstrates that the effect of increasing value of slip parameter 𝐾 is to shift the streamlines toward stretching boundary and thereby reducing thickness of the momentum boundary layer. Therefore, the effect of slip parameter 𝐾 is seen to decrease the boundary layer velocity while the temperature and concentration are increased with the increase of the slip parameter. The velocity curves show that the rate of transport decreases with the increasing distance (𝜂) normally to the sheet. In all cases, the velocity vanishes at some large distance from the sheet (at 𝜂=6).

The variations in velocity field, temperature distribution, and nanoparticle concentration profile for various values of 𝑀 are presented in Figures 5, 6, and 7. It is clear from these figures that the velocity decreases, whereas the temperature and concentration increase with the increase of the magnetic field parameter. This is due to the fact that the application of transverse magnetic field in an electrically conducting fluid produces a resistive force known as Lorentz force. This force has the tendency to slow down the motion of the fluid in the boundary layer and thus reducing the velocity. The additional work done in dragging the conducting nanofluid against the action of the magnetic field, 𝐵_𝑜, is manifested as thermal energy. This heats the conducting nanofluid and elevates temperatures. Thus, the presence of magnetic field decreases the momentum boundary layer thickness and increases the thermal boundary layer thickness. The warming of the boundary layer therefore also aids in nanoparticle diffusion which causes a rise in nanoparticle volume fraction, 𝜙.

Figures 8, 9, and 10 show the effect of viscoelastic parameter 𝛼 on the evolution of fluid motion and subsequently on the distribution of heat and mass across the sheet as time evolves. From this plot, it is evident that increasing values of viscoelastic parameter 𝛼 opposes the motion of the liquid close to the stretching sheet and assists the motion of the liquid far away from the stretching sheet. Increasing values of second-grade parameter enables the liquid to flow at a faster rate, due to which there is decline in the heat transfer. This is responsible for the increase in momentum boundary layer, whereas the thermal and concentration boundary layers reduce when the viscoelastic effects intensify.

The variation of nondimensional temperature profile 𝜃 with different values of heat source/sink parameter 𝑄 is shown in Figure 11. Increasing the heat source or sink parameter 𝑄 has tendency to increase the thermal state of the fluid. This increase in the fluid temperature causes more induced flow towards the plate through the thermal buoyancy effect. So, the thickness of the thermal boundary layer reduces for the increase of heat sink parameter, but it increases with heat source parameter.

In order to understand the influence of thermophoresis parameter Nt on heat and mass transfer, the local Nusselt and Sherwood numbers are plotted in Figures 12 and 13 for different values of slip parameter 𝐾. It is clear from these figures that the dimensionless heat transfer rates decrease with the increase in thermophoresis parameter Nt or slip parameter 𝐾, whereas the dimensionless mass transfer rates increase with the increase in thermophoresis parameter Nt. The numerical values of reduced Nusselt and Sherwood numbers, for different values of viscoelastic parameter 𝛼, against Nt are displayed in Figures 14 and 15. From the earlier graphical results, we have noticed that the thickness of the thermal and concentration boundary layers reduces when the viscoelastic effects intensify. This reduction is compensated with the increase in the rate of heat and mass transfer at the stretching surface. Also, the profiles obtained from these figures of reduced Nusselt and Sherwood numbers are compared with those reported by Noghrehabadi et al. [15] in case of 𝛼=0.

In the present study, the local rate of heat transfer (reduced Nusselt number Nur) and local rate of mass transfer at the sheet (reduced Sherwood number Shr), defined in (9), are the important characteristics. The numerical values of reduced Nusselt number and reduced Sherwood number are exhibited in Tables 1–3. Table 1 shows the excellent correlation between the current FEM computations and the earlier results of Khan and Pop [5], Wang [21], and Reddy Gorla and Sidawi [22], for reduced Nusselt number (−𝜃^′(0)) by neglecting the Brownian effect (Nb), slip effect (𝐾), magnetic effect (𝑀), and thermophoresis (Nt) for various values of Prandtl number (Pr) with step size, h = 0.0003. Variations of reduced Nusselt number Nur and reduced Sherwood number Shr with Nt are depicted in Table 2. Further, for various values of Nt and Nb, the results are compared with those reported by Khan and Pop [5] and Noghrehabadi et al. [15], and the comparison is found to be in good agreement for each value of Nt and Nb, which confirms that the present results are accurate. The numerical values of reduced Nusselt number and reduced Sherwood number are provided in Tables 3(a) and 3(b) for some values of Nt and Nb when Pr=Le=10,𝛼=0.5,𝑀=1.0,𝑄=0.1, and 𝐾=0.0. It is observed that reduced Nusselt number is a decreasing function, while reduced Sherwood number is an increasing function of Nt.

Table 1 

Comparison of results for reduced Nusselt number −𝜃^′(0).

Table 2 

Comparison of results for the reduced Nusselt number (−𝜃^′(0)) and reduced Sherwood number (−𝜙^′(0)) when Le=Pr=10 and 𝑀=𝑄=𝐾=𝛼=0.

Table 3 

(a) Variations of Nur with Nb and Nt when Pr=Le=10, 𝛼=0.5, 𝑀=1.0, 𝑄=0.1, and 𝐾=0.0. (b) Variations of Shr with Nb and Nt when Pr=Le=10, 𝛼=0.5, 𝑀=1.0, 𝑄=0.1, and 𝐾=0.0.

5. Conclusions

The problem of MHD boundary layer flow of a viscoelastic nanofluid past a stretching sheet has been solved numerically to exhibit the effect of partial slip (i.e., Navier’s condition) and heat source/sink on the fluid flow and heat transfer characteristics. By using a similarity transformation, the conservative equations for mass, momentum, energy, and nanoparticle concentration are transformed into the ordinary differential equations. We used the finite element method (FEM) for the numerical solution of these equations. The effects of different controlling parameters on the flow field and heat transfer characteristics are examined. The variation of reduced Nusselt number and Sherwood number with Nt for various values of 𝛼 and 𝐾 is presented in graphical form. The result can be summarized as follows.(1)With the increase in the second-grade parameter 𝛼, the velocity and the momentum boundary layer thickness increases; however, the temperature and nanoparticles concentration decrease.(2)Temperature increases and thermal boundary layer thickens when the strength of heat source/sink parameter 𝑄 increases.(3)There is a decrease in the velocity, but temperature and concentration are found to increase with an increase in velocity slip parameter 𝐾.(4)Magnetic field decelerates the flow and enhances temperatures and nanoparticle volume fraction (concentration) distributions in the boundary layer.(5)With the increase in the slip parameter 𝐾, heat transfer rate and mass transfer rate decrease.(6)By the increase of thermophoretic number Nt, the effect of velocity slip parameter 𝐾 on reduced Nusselt number Nur and reduced Sherwood number Shr increases and decreases, respectively.(7)The reduced Nusselt number and reduced Sherwood number both increases with the increase of viscoelastic parameter 𝛼.