Int. J. Appl. Comput. Math DOI 10.1007/s40819-016-0184-5 ORIGINAL PAPER

Finite Element Analysis of Viscoelastic Nanofluid Flow with Energy Dissipation and Internal Heat Source/Sink Effects P. Rana1 · R. Bhargava2 · O. Anwar Bég3 · A. Kadir3

© Springer India Pvt. Ltd. 2016

Abstract A numerical study is conducted of laminar viscoelastic nanofluid polymeric boundary layer stretching sheet flow. Viscous dissipation, surface transpiration (suction/injection), internal heat generation/absorption and work done due to deformation are incorporated using a second grade viscoelastic non-Newtonian nanofluid with non-isothermal associated boundary conditions. The nonlinear boundary value problem is solved using a higher order finite element method. The influence of viscoelasticity parameter, Brownian motion parameter, thermophoresis parameter, Eckert number, Lewis number, Prandtl number, internal heat generation and also wall suction on thermofluid characteristics is evaluated in detail. Validation with earlier non-dissipative studies is also included. The hp-finite element method achieves the desired accuracy at p = 8 with comparatively less CPU cost per iteration (with less degrees of freedom, DOF) as compared to lower order finite element methods. The simulations have shown that greater polymer fluid viscoelasticity (k1 ) accelerates the flow. A rise in Brownian motion parameter (Nb) and thermophoresis parameter (Nt) elevates temperatures and reduce the heat transfer rates (local Nusselt number function). Increasing Eckert number increases temperatures whereas increasing Prandtl number (Pr) strongly lowers temperatures. Increasing internal heat generation (Q > 0) elevates temperatures and reduces the heat transfer rate (local Nusselt number function) whereas heat absorption (Q < 0) generates the converse effect. Increasing suction ( f w > 0) reduces velocities and temperatures but elevates enhances mass transfer rates (local Sherwood number function), whereas increasing injection ( f w < 0) accelerate the flow, increases temperatures and depresses wall mass transfer rates. The study finds applications in rheological nano-bio-polymer manufacturing.

B

O. Anwar Bég [email protected]; [email protected]

1

Department of Mathematics, Jaypee Institute of Information Technology, Noida, India

2

Department of Mathematics, Indian Institute of Technology, Roorkee 247667, India

3

Spray Research Group, School of Computing, Science and Engineering, Newton Bldg, The Crescent, University of Salford, Manchester M54WT, England, UK

123

Int. J. Appl. Comput. Math

Keywords Viscoelastic nanofluid · Eckert number · Lewis number · Stretching · Brownian motion · Thermophoresis · hp-FEM · Lewis number · Nano-polymer materials processing

List of symbols Roman uw A, B, E Q C Cw C∞ Nt (x, y) Tw T∞ T Pr qm qw DB DT uw f (η) g(η) Nb Le k1 N ux A1 , A2 Ec Sh x Cf u, υ fw m

Sheet velocity (m/s) Constants Internal heat source/sink Nanoparticle volume fraction Nanoparticle volume fraction Ambient nanoparticle volume fraction Thermophoresis parameter Cartesian coordinates (m) Temperature at the sheet (K) Ambient temperature attained (K) Temperature on the sheet (K) Prandtl number Wall mass flux (kg/s) Wall heat flux (W/m2 ) Brownian diffusion coefficient (m2 /s) Thermophoretic diffusion coefficient (m2 /s) Velocity of stretching sheet (m/s) Dimensionless stream function Gravitational acceleration (m/s2 ) Brownian motion parameter Lewis number Viscoelastic parameter Nusselt number Rivlin–Ericksen tensors in the constitutive relation (N/m2 ) Eckert number Sherwood number Skin friction Velocity components along x–y axes (m/s) Suction/injection parameter Power-law parameter

Greek symbols  τ (ρc) f φ (η) η θ (η)

123

Stress tensor (N/m2 ) Parameter defined by ε (ρc) p / (ρc) f Heat capacity of the fluid (J/kg3 K) Rescaled nanoparticle volume fraction Similarity variable Dimensionless temperature

Int. J. Appl. Comput. Math

(ρc) p ρf β ρp ψ ν αm α1 , α2 β1 , β2 , β3

Effective heat capacity of the nanoparticle material (J/kg3 K) Fluid density (kg/m3 ) Volumetric expansion coefficient of the fluid (l/K) Nanoparticle mass density (kg/m3 ) Stream function Fluid kinematic viscosity (m2 /s) Thermal diffusivity (m2 /s) Material moduli (N/m2 ) Higher order viscosities (m2 /s)

Subscripts w ∞

Condition on the sheet (wall) Condition far away from the sheet (free stream)

Introduction In recent years, an important trend in biopolymer manufacture has been the facility of patterning functional materials at different length scales [1]. This characteristic allows the precision fabrication of functional biopolymers for many diverse applications including cell biology, tissue engineering, bio-optics (contact lenses, ophthalmic agents etc) and suspension fluids for medical transplants e.g. biological hydrogels [2]. Many of these new biopolymers exhibit beneficial rheological properties enabling enhanced performance in delivery of alginates in for example the accelerated treatment of vascular hemorrhages, arteriovenous malformations via injection in medical micro-catheters for endovascular embolization [3]. The necessity to enhance mass transfer (oxygen) in microbial biopolymers [4] has also attracted the science of nanotechnology to the optimization of such materials. As a result nano-bio-polymeric fluid dynamics has emerged as an exciting new research area in medical engineering. This domain combines the properties of functional biopolymers and nanofluids to achieve better and more adaptive agents for treatments via enhanced heat and mass transfer features. Nanofluids [5,6] describe a solid–liquid mixture which consists of a fluid suspension containing ultra-fine particles termed nanoparticles. The nanoparticles Al2 O3 , CuO, TiO2 , ZnO and SiO2 are increasingly being employed in biomedical systems. Typical thermal conductivity enhancements for bio-nanofluids [7] are in the range of 15–40 % over the base fluid and heat transfer coefficients enhancements have been found up to 40 %. Pak and Cho [8] conducted comparative experimental investigations on turbulent friction and heat transfer of nanofluids with alumina and titanium oxide nanoparticles in a circular pipe. They tested alumina and titanium oxide nanoparticles with mean diameters of 13 and 27 nm respectively, in water. They found that inclusion of a 10 % volume fraction of alumina in water increased the viscosity of the fluid 200 times and inclusion of the same volume fraction of titanium oxide produced a viscosity that was 3 times greater than water. They also noticed the heat transfer coefficient increased by 45 % at 1.34 % volume fraction of alumina to 75 % at 2.78 % volume fraction of alumina and is consistently higher than titanium oxide nanofluid. Increases in thermal conductivity of this magnitude in nanofluids cannot be solely attributed to the higher thermal conductivity of the added nanoparticles and therefore other mechanisms and factors must contribute. These include particle agglomeration [9–11], volume fraction [9], Brownian motion [12,13], thermophoresis, nanoparticle size [13], particle shape/surface area [2,14],

123

Int. J. Appl. Comput. Math

liquid layering on the nanoparticle-liquid interface [15], temperature [13,16] and reduction of thermal boundary layer thickness. The literature on the thermal conductivity and viscosity of nanofluids has been reviewed by Eastman et al. [17], Wang and Mujumdar [18] and Trisaksri and Wongwises [19]. In addition, a succinct review on applications and challenges of nanofluids has also been provided by Wen et al. [20] and Saidur et al. [21]. Buongiorno [16] identified multiple mechanisms in the convective transport in nanofluids using a two-phase non-homogenous model including inertia, Brownian diffusion, thermophoresis, diffusiophoresis, the Magnus effect, fluid drainage and gravity. Of all of these mechanisms, only Brownian diffusion and thermophoresis were found to be important in the absence of turbulence effects. He also suggested that the boundary layer has different properties owing to the effect of temperature and thermophoresis. Taking Brownian motion and thermophoresis into account, he developed a correlation for the Nusselt number which was compared to data from Pak and Cho [8] and which correlated best with the latter [8] experimental data. Recently, the Buongiorno [16] model has been used by Kuznetsov and Nield [22] to study the natural convection flow of nanofluid over a vertical plate and their similarity analysis identified four parameters governing the transport process. Rana et al. [23] investigated the mixed convection problem along an inclined plate in the porous medium. In the alternative approach, the heat transfer analysis with non-uniform heating along a vertical plate has been studied by Rana and Bhargava [24]. Prasad et al. [25] studied micropolar nanofluid convection from a cylinder using a finite difference scheme. Khan and Pop [26] used the Kuznetsov–Nield model to study the boundary layer flow of a nanofluid past a stretching sheet with a constant surface temperature. Subsequently several authors have solved the problem of nanofluid flow from a stretching sheet with different stretching and boundary conditions and representative articles in this regard include Rana and Bhargava [27], Uddin et al. [28,29], Nedeem and Lee [30], Kandasamy et al. [31], Bachok et al. [32] and Rana et al. [33]. As mentioned earlier non-Newtonian (rheological) behaviour of nanofluids has been highlighted in nano-polymeric manufacturing processes by many researchers including Chen et al. [34], Chen et al. [35] and Gallego et al. [36]. Recently, Khan and Gorla [37] investigated heat and mass transfer in non-Newtonian nanofluid over a non-isothermal stretching wall. Numerous applications of viscoelastic nanocomposite fabrication techniques [3,4] have led to renewed interest among researchers to investigate viscoelastic boundary layer flow over a stretching sheet. Although numerous studies of viscoelastic stretching flows have been communicated by, for example, Rajagopal et al. [38,39], Dandapat and Gupta [40] and Rao [41], these do not study nanotechnological materials. To improve our understanding of the manufacture of bio-nano-polymers via e.g. sheet extrusion from a dye, it is important to study heat transfer (especially cooling rates) which strongly influences the constitution and quality of manufactured products in biomaterials processing. Viscoelastic fluid flow generates heat by means of viscous dissipation and work done due to deformation. There is another important aspect, which should also be taken into the account in a situation when there would be a temperature-dependent heat source/sink present in the boundary layer region. A wide variety of problems with heat and fluid flow over a stretching sheet have been studied with viscoelastic fluids and with different thermal boundary conditions (prescribed surface temperature, PST and prescribed heat flux, PHF) and power-law variation of the stretching velocity. A representative sample of recent literature on computational simulations of viscoelastic flows based on Reiner–Rivlin differential models is provided in [42–46]. To the authors’ knowledge very few studies have thus far been communicated with regard to boundary layer flow and heat transfer of a viscoelastic polymeric nanofluid (Ethylene glycol and polymer based nanofluids) extruding from a stretching sheet with energy dissipation. This problem is very relevant to modern nano-polymeric fluid manufacture processes

123

Int. J. Appl. Comput. Math

in the biotechnology industries wherein materials can be synthesized for specific medical applications including sterile coating, anti-bacterial buffers etc. [1–4]. However some recent efforts to simulate viscoelastic flows have been communicated with alternative viscoelastic constitutive equations to the second grade differential model. For example Krishnamurthy et al. [47] employed the Williamson viscoelastic model to study reactive phase change heat transfer in nanofluid dynamics of porous media. Hussain et al. [48] employed the Jefferys viscoelastic model for magnetized Sakiadis nanofluid flow with exponential stretching and thermal radiative effects. Khan et al. [49] deployed an Oldroyd-B viscoelastic model to study rheological effects on heat and mass transfer in 3-dimensional nanofluid boundary layer flow. Akbar et al. [50] applied the elegant Eyring–Powell rheological model for magnetohydrodynamic stretching sheet flow. Mehmood et al. [51] also used the Jefferys model for oblique stagnation flow and heat transfer. Haq et al. [52] used the Eringen microplar model to study nanofluid stagnation point flow with free convection and radiative effects. All these studies demonstrated the significant influence of viscoelasticity in modifying velocity and also heat and mass transfer distributions. In the present paper, we focus on refining the simulations for bio-rheological polymer materials by incorporating a proper sign for the normal stress modulus (i.e., α ≥ 0) as described in “Constitutive Relations and Bio-Nano-Polymer Mathematical Model” section. A numerical solution is developed for the nonlinear boundary value problem derived in the present article. Babuska and Guo [53] presented the basic theory and applications of h, p and h- p versions of the finite element methods. Khomami et al. [54] has presented a comparative study of higher and lower order finite element techniques for computation of viscoelastic flows. It has been demonstrated that the hp-finite element method gives rise to an exponential convergence rate toward the exact solution, while all the lower order schemes considered exhibit a linear convergence rate. Thus, we employ an extensively validated, highly efficient, hp-Galerkin finite element method to obtain numerical solutions for the present problem. This paper runs as follows. In “Constitutive Relations and Bio-Nano-Polymer Mathematical Model” section, we consider the mathematical analysis of the viscoelastic polymeric nanofluid flow and heat transfer with energy dissipation effects (viscous heating and work done due to deformation) cited above are included in the energy equation with prescribed surface temperature (PST case) boundary heating. A summary of the hp-finite element numerical technique is presented in “Numerical Solution” section. “Results and Discussion” section presents graphical solutions and a discussion of the influence of the non-dimensional parameters including Brownian motion parameter, thermophoresis parameter, viscoelastic parameter, Prandtl number, Eckert number, the surface suction/injection parameter, and internal heat source/sink parameter, on the flow characteristics. Finally the conclusions follow in “Conclusions” section.

Constitutive Relations and Bio-Nano-Polymer Mathematical Model Let us recall the constitutive equation for an incompressible fluid of grade n (based on the postulate of gradually fading memory) given by Coleman and Noll [55]:  (t) = −P I +

n 

Sj

(1)

j=1

For n = 3, the first three tensors S j are given by: S1 = μA1 .

123

Int. J. Appl. Comput. Math

S2 = α1 A2 + α2 A21 .

  S3 = β1 A3 + β2 (A2 A1 + A1 A2 ) + β3 trA21 A1 .

(2)

where  is the stress tensor, −PI designates the indeterminate part of the stress, μ is the viscosity, α1 , α2 are the normal stress moduli, and β1 , β2 , β3 are the higher order viscosities. The Rivlin–Ericksen tensors An are defined by the recursion relation: A1 = (grad V ) + (grad V )T . d An = dt An−1 + An−1 . (grad V ) + (grad V )T .An−1 , n = 2, 3, . . .

(3)

where V denotes the velocity field, grad is the gradient operator and d/dt is the material time derivative. Thus, for the particular of a second-grade fluid, we have:  = − p I + μA1 + α1 A2 + α2 A21

(4)

where  is the Cauchy stress tensor, p is the pressure, μ is the viscosity, α1 and α2 are two normal stress moduli with α1 < 0, and A1 and A2 are the first two Rivlin–Ericken tensors defined by: A1 = (grad V ) + (grad V )T . (5) A2 = ddtA1 + A1 . (grad V ) + (grad V )T .A1 The model (4) displays normal stress differences in shear flow and is an approximation to simple fluid in the sense of retardation. This model is applicable to some dilute polymer solutions and is valid at low rates of shear. Dunn and Fosdick [56] have shown that, for the fluid modelled by Eq. (4), to be compatible with thermodynamics and to satisfy the Clausius– Duhem inequality for all motions, and the assumption that the specific Helmholtz free energy of the fluid takes its minimum values in equilibrium, the material moduli must satisfy: μ ≥ 0, α ≥ 0, α1 + α2  = 0.

(6)

However for many of the non-Newtonian fluids of rheological interest, the experimental results for α1 and α2 do not satisfy the restrictions (6). By using data reduction from experiments, Fosdick and Rajagopal [57] have shown that in the case of a second-order fluid the material moduli, μ, α1 and α2 should satisfy the following relations: μ ≥ 0, α1 ≤ 0, α1 + α2  = 0.

(7)

They also found that the fluids modeled by Eq. (4) with the relationship (7) exhibit some anomalous behaviors. A critical review on this controversial issue can be found in the work of Dunn and Rajagopal [58]. Generally, in the literature the fluid obeying Eq. (4) with α1 < 0 is termed as a second-order fluid and with α1 > 0 is termed as second-grade fluid. When α1 = 0, α2 = 0 and μ > 0 Eq. (4) reduces to the well-known constitutive relation for an incompressible Newtonian fluid. Equation (4) is used in the present simulation. We consider steady, incompressible, laminar, two-dimensional boundary layer flow of a viscoelastic polymeric nanofluid past a flat sheet coinciding with the plane y = 0 with the flow being confined to y > 0. The flow is generated, due to non-linear stretching of the sheet, caused by the simultaneous application of two equal and opposite forces along the x-axis. Keeping the origin fixed, the sheet is then stretched with a velocity u w = E x where E is a constant, and x is the coordinate measured along the stretching surface, varying linearly with the distance from the slit. A schematic representation of the physical model and coordinates system is depicted in Fig. 1. The pressure gradient and external forces are neglected. The stretching surface is maintained at constant temperature and concentration, Tw and Cw respectively, and

123

Int. J. Appl. Comput. Math

Fig. 1 Physical model and co-ordinate system

these values are assumed to be greater than the ambient temperature and concentration, T∞ and C∞ , respectively. The governing equations for conservation of mass, momentum, thermal energy and nanoparticle species diffusion, for a second order Reiner–Rivlin viscoelastic nanofluid can be written in Cartesian coordinates, x, y as: ∂u ∂v + =0 (8) ∂x ∂y      2  ∂ ∂ 3υ ∂u ∂u ∂ 2 υ ∂T ∂ 2u ∂ u + υ ρ u (9) +υ = μ 2 + α1 u 2 + ∂x ∂y ∂y ∂x ∂y ∂ y ∂ y2 ∂ y3

    ∂T 2 ∂C ∂ T ∂T ∂T ∂2T ν ∂u 2 u +υ = α 2 + τ DB . + (DT /T∞ ) + ∂x ∂y ∂y ∂y ∂y ∂y c f ∂y   α1 ∂u ∂ ∂u ∂u q (T − T∞ ) + (10) +ν u + ∂x ∂y (ρc) f ∂ y ∂ y (ρc) f u

∂2T ∂ 2C ∂C ∂C +υ = D B 2 + (DT /T∞ ) 2 ∂x ∂y ∂y ∂y

where αm =

(ρc) p km , τ= (ρc) f (ρc) f

(11)

(12)

subject to the boundary conditions  x 2 , υ = υw , u w = E x, T = Tw (x) = T∞ + A l  x 2 C = Cw (x) = C∞ + B at y = 0 l u = υ = 0, du/dy = 0, T = T∞ , C = C∞ as y → ∞

(13a) (13b)

Here u and v are the velocity components along the axes x and y, respectively, α1 is the modulus of the viscoelastic fluid, ρ f is the density of the base fluid, αm is the thermal diffusivity, υ is the kinematic viscosity, E is a positive constant, D B is the Brownian diffusion coefficient, DT is the thermophoretic diffusion coefficient and τ = (ρc) p / (ρc) f is the ratio

123

Int. J. Appl. Comput. Math

between the effective heat capacity of the nanoparticle material and heat capacity of the fluid, c is the volumetric volume expansion coefficient and ρ p is the density of the particles. Equations (9)–(13) are a new formulation and extend the Newtonian model of Rana and Bhargava [27] to a second order model, by incorporating new terms from the Reiner–Rivlin model and also Eckert heating, heat generation/absorption and work deformation terms. Proceeding with the analysis, we introduce the following dimensionless variables: η = (b/ν)1/2 y, u = bx f  (η) , ψ = (bν)1/2 y  x 2  x 2 θ (η) , C − C∞ = B φ (η) (14) T − T∞ = A l l where the stream function ψ is defined in the usually way as u = ∂ψ/∂ y and υ = −∂ψ/∂ x. In seeking a similarity solution based on the transformations in Eq. (14), we have taken into account that the pressure in the outer (inviscid) flow is p = p0 (constant). The governing Eqs. (8)–(11) then reduce to:

  2 (15) f  + f f  − f 2 + k1 2 f  f  − f  − f f I V = 0  2 1  θ + f θ  + N bθ  φ  + N t θ  − 2 f  θ + Qθ Pr

    2 +Ec f  + k1 f  f  f  − f f  = 0

(16)

  N t  φ  + Le f φ  − 2 f  φ + θ =0 Nb The transformed boundary conditions are η = 0, f = f w , η → ∞,

(17)

f  = 1, θ = 1, φ = 1

(18a)

f  = 0, θ = 0, φ = 0

(18b)

where () denotes differentiation with respect to η and the key dimensionless thermo-physical parameters are defined by: Pr =

ν , α

Ec =

E2 L2 , Ac p

Le =

ν α1 E , , k1 = DB ρν

Nb =

f w = −υw / (Eν)1/2 ,

(ρc) p D B (Cw − C∞ ) , (ρc) f ν

Nt =

Q=

q E (ρc) f

(ρc) p DT (Tw − T∞ ) (ρc) f νT∞

(19)

Here Pr, k1 , Le, N b, N t, Ec, f w and Q denote the Prandtl number, viscoelastic parameter, Lewis number, Brownian motion parameter, thermophoresis parameter, Eckert number, surface suction/injection parameter and internal heat source/sink parameter, respectively. It is important to note that this boundary value problem reduces to the problem of flow and heat and mass transfer due to a stretching surface in a viscoelastic fluid when N b, N t are zero in Eqs. (16) and (17). The presence of viscoelastic terms in the momentum equation (15) raises the order of this equation to one above that of the Navier–Stokes equations. Well-posedness of the problem can be achieved via a number of strategies. To give a general review of the past work on the existence of solutions of the primitive Eqs. (8) and (9), Troy et al. [59] and Mcleod and Rajagopal [60] obtained a unique solution. In fact, the second condition in Eq. (13b) is the property of the boundary layer in the asymptotic region. Chang [61] has claimed that the solution of the problem is not necessarily unique without this condition. In the absence of the viscoelastic parameter (i.e. k1 = 0), Eq. (9) reduces to a third-order ordinary differential equation for which these four conditions are also applicable. There is further an analytical solution in the absence of slip and viscoelasticity, which reads:

123

Int. J. Appl. Comput. Math

f (η) =

1 − e−mη 1 + f w with m = m 2

 fw +



 4 + 4 f w2

(20)

Among all these, solutions of the form proposed by Troy et al. [59] are the realistic ones as we can recover the boundary layer approximation of Navier–Stokes solution only in the limiting case of k1 = 0. It is worth mentioning that Eq. (9) with the boundary conditions Eqs. (13a, 13b) has an exact solution as given by: f (η) =

1 − e−mη + fw , m

(21)

where m is a real positive root of the cubic algebraic equation: k1 f w m 3 + (k1 + 1) m 2 − f w m − 1 = 0

(22)

The velocity profile is determined from Eq. (21) to be: f  (η) = e−mη The skin friction coefficient C f can further be determined from:  2 ∂ u ∂u ∂υ u + α − 2 μ ∂u 1 ∂y ∂ x∂ y ∂y ∂y Cf = ρu 2w /2

(23)

(24)

Noting that the skin friction parameter is − f  (0) = r , Eq. (24) can be further simplified to:  b x (25) C f = −2m (1 + 3k1 ) υ For the linearly stretching boundary layer problem, the exact solution for f is f (η) = 1−e−η and this exact solution is unique. Important heat and mass transfer quantities of practical interest for the present flow problem are the local Nusselt number and the local Sherwood number which are defined, respectively as: xqw xqm N ux = , Sh x = (26) k (Tw − T∞ ) D B (Cw − C∞ ) where qw and qm are heat flux and mass flux at the surface (plate), respectively. The dimensionless heat and mass transfer rates can also be shown to take the form given in following expressions:   υ υ  , (27) N u = −θ Sh x = −φ  (0) (0) x 2 ax ax 2 The set of ordinary differential equations defined by Eqs. (15)–(17) are highly non-linear and cannot be solved analytically. The hp-finite element method [53] has therefore been implemented to solve this highly coupled two-point boundary value problem to determine the velocity, temperature and nano-particle concentration distributions.

Numerical Solution The Finite Element Method Finite element method (FEM) was basically developed in reference to aircraft structural mechanics problems and has evolved over a number of decades to become the dominant

123

Int. J. Appl. Comput. Math

computational analysis tool for solving the linear and non-linear ordinary differential, partial differential and integral equations. The finite element method provides superior versatility to other numerical methods include finite differences and is generally very stable with excellent convergence characteristics.

Finite-Element Discretization The whole domain is divided into a finite number of sub-domains, designated as the discretization of the domain. Each sub-domain is called an element. The collection of elements comprises the finite-element mesh.

Generation of the Element Equations From the mesh, a typical element is isolated and the variational formulation of the given problem over the typical element is constructed. An approximate solution of the variational problem is assumed and the element equations are generated by substituting this solution in the above system. The element matrix, which is called stiffness matrix, is constructed by using the element interpolation functions.

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

Imposition of Boundary Conditions The essential and natural boundary conditions are imposed on the assembled equations.

Solution of Assembled Equations The assembled equations so obtained can be solved by any of numerical technique including the Gauss elimination method, LU Decomposition method, Householder’s technique, Choleski decomposition etc.

hp-Finite Element Method hp-FEM is a general version of the finite element method (FEM), based on piecewisepolynomial approximations that employ elements of variable size (h) and polynomial degree ( p). The origins of hp-FEM date back to the pioneering work of Babuska et al. [53] who discovered that the finite element method converges exponentially faster when the mesh is refined using a suitable combination of h−refinements (dividing elements into smaller ones) and p-refinements (increasing their polynomial degree). The exponential convergence makes the method a very attractive choice compared to most other finite element methods which only converge at an algebraic rate. An excellent demonstration of the exponential convergence rate of hp-FEM for viscoelastic fluid flows has been provided by Khomami et al. [54]. For

123

Int. J. Appl. Comput. Math

the solution of system of simultaneous, coupled, nonlinear systems of ordinary differential equations as given in (15–17), with the boundary conditions (18), we first assume: ∂f =h ∂η

(28)

The system of Eqs. (8–10) then reduces to

 2  ∂h 2n ∂ 3h ∂ 2h 2 − f 3 =0 h + k1 2h 2 − n+1 ∂η ∂η ∂η   ∂θ ∂φ ∂θ 2 ∂θ 1 ∂ 2θ + N b + N t + f − 2hθ + Qθ Pr ∂η2 ∂η ∂η ∂η ∂η  

  ∂h 2 ∂h ∂ 2h ∂h + Ec + k1 − f 2 h =0 ∂η ∂η ∂η ∂η   ∂φ ∂f N t ∂ 2θ ∂ 2φ + Le f =0 + 2 φ + ∂η2 ∂η ∂η N b ∂η2 

∂ 2h ∂h + f − ∂η ∂η

(29)

(30) (31)

and the corresponding boundary conditions now become; f = fw ,

η = 0, η → ∞,

f  = 1, θ = 1, φ = 1



f = 0, θ = 0, φ = 0

(32a) (32b)

Variational Formulation The variational form associated with Eqs. (28)–(31) over a typical linear element, e = (ηe , ηe+1 ), is given by: ηe+1



w1 ηe ηe+1



w2 ηe ηe+1

 ∂f − h dη = 0 ∂η ∂ 2h ∂h + f − ∂η ∂η



(33)

  2  ∂h 2n ∂ 3h ∂ 2h 2 − f 3 dη = 0 h + k1 2h 2 − n+1 ∂η ∂η ∂η (34)





2

1 ∂ 2θ ∂θ ∂φ ∂θ ∂θ + Nb + Nt + f − 2hθ + Qθ Pr ∂η2 ∂η ∂η ∂η ∂η ηe      ∂h 2 ∂h ∂h ∂ 2h h + k1 + Ec − f 2 dη = 0 ∂η ∂η ∂η ∂η w3

ηe+1

w4 ηe



∂ 2φ + Le ∂η

 f

  ∂φ ∂f N t ∂ 2θ dη = 0 +2 φ + ∂η ∂η N b ∂η2

(35)

(36)

where w1 , w2 , w3 and w4 are arbitrary test functions and may be viewed as the variation in f, h, θ and φ respectively.

123

Int. J. Appl. Comput. Math

Finite Element Formulation The finite element model may be obtained from above equations by substituting finite element approximations of the form; f =

p 

f jψj, h =

j=1

p 

h jψj, θ =

j=1

p 

θjψj, φ =

p 

j=1

φjψj

(37)

j=1

with w1 = w2 = w3 = w4 = ψi , (i = 1, 2, . . . , p)

(38)

In our computations, the shape functions for a typical element () are:

.... x p

x1 x2 x3 x4 In global coordinates:

  (x − x1 ) (x − x2 ) (x − x4 ) . . . x − x p   i = 1, . . . , p ψi = (xi − x1 ) (xi − x2 ) (xi − x4 ) . . . xi − x p

(39)

In local coordinates: For p = 2 (linear element)

he

ηe ψ1e =

η e+1

(ηe+1 − η) (η − ηe ) , ψ2e = , ηe ≤ η ≤ ηe+1 (ηe+1 − ηe ) (ηe+1 − ηe )

p = 3 (Quadratic element)

he

ηe ψ1e = ψ3e =

(40)

(ηe+1 +ηe −2η)(ηe+1 −η) , (ηe+1 −ηe )2 (ηe+1 +ηe −2η)(η−ηe ) − , (ηe+1 −ηe )2

η e+1 ψ2e =

4(η−ηe )(ηe+1 −η) , (ηe+1 −ηe )2

ηe ≤ η ≤ ηe+1

(41)

p = 4 (Cubic element) he

ηe η e+1 ψ1e = ψ3e =

(ηe+1 +2ηe −3η)(2ηe+1 +ηe −3η)(η−ηe+1 ) e+1 −η) , ψ2e = 9(η−ηe )(2ηe+1 +ηe −3η)(η , 2(ηe+1 −ηe )3 2(ηe+1 −ηe )3 9(η−ηe )(ηe+1 +2ηe −3η)(ηe+! −η) (η−ηe )(ηe+1 +2ηe −3η)(ηe+1 −η) e , ψ4 = , 2(ηe+1 −ηe )3 2(ηe+1 −ηe )3

ηe ≤ η ≤ ηe+1 and so on.

123

(42)

Int. J. Appl. Comput. Math

The finite element model of the equations thus formed is given by; ⎡  11   12   13   14  ⎤ ⎡ ⎤ ⎡ b 1  ⎤ K K K K { } f ⎢ b 2  ⎥ ⎢  K 21   K 22   K 23   K 24  ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ {g} ⎥ ⎥        = ⎢ ⎥ ⎢ 31 ⎥ K 32 K 33 K 34 ⎦ ⎣ {θ } ⎦ ⎣ b3 ⎦ ⎣ K  41   42   43   44    {φ} K K K K b4

(43)

    where K mn and bm (m, n = 1, 2, 3, 4) are defined as: ηe+1

K i11 j

= ηe

ηe+! K i21 j = −

ψi ηe ηe+1

K i22 j

ηe+1

∂ψ j ψi dη, K i12 j =− ∂η

=− ηe

14 ψi ψ j dη, K i13 j = K i j = 0, ηe

∂h 24 ψ j dη, K i23 j = K i j = 0, ∂η

∂ψi ∂ψ j ∂η − ∂η ∂η

⎡

η! e+1

∂ψi ∂η

ηe+1

¯ j dη ψi hψ

ηe η! e+1

∂ψ h¯ ∂ηj

" ∂ψi ∂η

⎤

¯ ψ j ∂∂ηh dη

dη + ⎥ ⎢ −2 ⎥ ⎢ ηe + k1 ⎢ ηe+1ηe ⎥ η η e+1 e+1 2ψ ! ∂ψi ∂ 2 ψ ! ! ⎦ ⎣ ¯ ∂ ∂ψ ∂f j j ∂h ¯ + ψi ∂η ∂η2 dη − ψi ∂η ∂η dη ∂η f ∂η2 dη + ηe ηe+1

K i31 j = −Pr Eck1

ψ

ηe ηe+1

ψi

+ Pr Ec ηe ηe+1

K i33 j = − ηe

ηe

∂h ∂ 2 h ψ j dη, K i32 j = Pr Eck1 ∂η ∂η2

ηe ηe+1

ψi ηe

∂h ∂h ψ j dη ∂η ∂η

∂h ∂ψ j ∂η ∂η

∂ψi ∂ψ j dη + Nr Pr ∂η ∂η

ηe+1

ψi θ  ηe

∂ψ j dη + Pr ∂η

ηe+1

ψi h¯

ηe

∂ψ j dη ∂η

ηe+1

¯ j ψi hψ

− 2Pr ηe ηe+1

ψi θ 

K i34 j = N b Pr ηe

∂ψ j dη, ∂η ηe+1

42 43 K i41 j = K i j = 0, K i j = −N t ηe ηe+1

K i44 j

= −N b ηe

∂ψi ∂ψ j dη, ∂η ∂η

∂ψi ∂ψ j dη + LeN b ∂η ∂η

ηe+1 ηe

∂ψi dη − 2N bLe ∂η

ηe+1

¯ j dη, ψi hψ

ηe

123

Int. J. Appl. Comput. Math Table 1 Calculation of Nusselt number and Sherwood number when Nb = 0.3, Nt = 0.3, Pr = 10, Le = 10, k1 = 1.0, Q = 0.5, Ec = 0.0, f w = 0 E

p

| − θ  (0) |

DOF

η∞ = 4

| − φ  (0) | η∞ = 6

η∞ = 8

η∞ = 4

Total CPU time (s) η∞ = 6

η∞ = 8

η∞ = 6

400

2

3204

0.9754

0.9782

0.9792

4.9215

4.8999

4.8867

36.61

1000

2

8004

0.9720

0.9730

0.9738

4.9545

4.9492

4.9402

96.08

2000

2

16,004

0.9708

0.9713

0.9718

4.9606

4.9658

4.9628

207.25

4000

2

32,004

0.9703

0.9707

0.9711

4.9771

4.9741

4.9733

455.16

8000

2

64,004

0.9702

0.9705

0.9707

4.9799

4.9782

4.9772

1229.52

10,000

2

80,004

0.9700

0.9702

0.9703

4.9804

4.9797

4.9789

1642.70

20,000

2

160,004

0.9698

0.9699

0.9699

4.9801

4.9793

4.9790

5833.80

500

4

10,004

0.9719

0.9732

0.9744

4.9546

4.9400

4.9362

92.94

500

6

12,004

0.9711

0.9715

0.9719

4.9778

4.9603

4.9588

231.33

500

8

16,004

0.9700

0.9701

0.9701

4.9802

4.9793

4.9789

1128.11

500

10

20,004

0.9699

0.9700

0.9700

4.9799

4.9792

4.9789

3811.23

Bold values corresponds to mesh-independent solution E, Number of elements; p, degree of polynomial; DOF, degrees of freedom

    ηe+1 df dθ ηe+1 ∂ 2h ∂h bi1 = 0, bi2 = − ψi , bi3 = − ψi , + k1 ψi d h¯ − f¯ 2 dη ∂η ∂η dη ηe ηe ηe+1   dθ dφ + (44) bi4 = − ψi dη dη ηe where h¯ =

p  i=1

h i ψi , h  =

p 

hi

i=1

 ∂ψi  ∂ψi ∂ψi φi θ¯ , θ = , φ = ∂η ∂η ∂η p

p

i=1

i=1

(45)

Validation of Numerical Solution Validation of the present numerical solutions is demonstrated in two ways. Firstly, an extensive mesh testing procedure for the h-type and hp-type Galerkin schemes has been conducted to ensure a grid-independent solution for the given boundary value problem, as documented in Table 1. It is observed that in the same domain by increasing the polynomial degree of approximation, one can achieve the desired accuracy with less DOF (Table 1). Also, arbitrary values of the thermophysical parameters are selected to verify the results and very little variation observed in the computations. The total CPU times on a Dell T5500 system have also been included in Table 1 and it is apparent that the desired accuracy is achieved for both heat and mass transfer rates, with p = 8 with an optimal time of 1128.11 s. Thus, for the present study a polynomial degree of approximation, p = 8 and number of elements, E = 500 have been adopted. Secondly, In order to verify the accuracy of the numerical solutions, the validity of the present numerical code has been benchmarked for the special case of Newtonian flow in the absence of viscous heating, heat source/sink, wall suction and vanishing thermophoresis and Brownian motion effects with constant surface temperature of the sheet i.e. with k1 = 0, Ec = 0, Q = 0, f w = 0, N b = N t = 10−5 and CST. This special

123

Int. J. Appl. Comput. Math Table 2 Comparison of results for the reduced Nusselt number, −θ  (0) with k1 = 0, Ec = 0, Q = 0, f w = 0, N b = N t = 10−5 and CST Pr

Wang [62]

Gorla and Sidawi [63]

Khan and Pop [26]

Present results

0.07

0.0656

0.0656

0.0663

0.0655

0.20

0.1691

0.1691

0.1691

0.1691

0.70

0.4539

0.5349

0.4539

0.4539

2.00

0.9114

0.9114

0.9113

0.9113

7.00

1.8954

1.8905

1.8954

1.8953

20.00

3.3539

3.3539

3.3539

3.3539

70.00

6.4622

6.4622

6.4621

6.4621

Table 3 Comparison of results for the reduced Nusselt number, −θ  (0) with k1 = 0, Ec = 0, Q = 0, f w = 0, Pr = Le = 10 and CST

Nb

Nt

Nur [26]

Shr [26]

Nur present

Shr present

0.1

0.1

0.9524

2.1294

0.9524

2.1294

0.2

0.1

0.5056

2.3819

0.5056

2.3819

0.3

0.1

0.2522

2.4100

0.2521

2.4101

0.4

0.1

0.1194

2.3997

0.1194

2.3999

0.5

0.1

0.0543

2.3836

0.0541

2.3836

0.1

0.2

0.6932

2.2740

0.6932

2.2740

0.1

0.3

0.5201

2.5286

0.5201

2.5286

0.1

0.4

0.4026

2.7952

0.4026

2.7952

0.1

0.5

0.3211

3.0351

0.3210

3.0352

case was studied earlier by Wang [62], Gorla and Sidawi [63] and Khan and Pop [26] and inspection of Table 2 shows excellent correlation of local Nusselt number (−θ  (0)) with CST as computed with hp-FEM and the other published computations, for different values of Pr . In Table 3, the hp-FEM results are further compared with Khan and Pop [26] for different combinations of Nb and Nt. The results are also validated with earlier computations of Nataraja et al. [64], Mushtaq et al. [65] and Chen [45] for second-grade viscoelastic fluid flow with PST keeping Ec = 0, Q = 0, k1 = 0, f w = 0, N b = N t = 10−5 and no work due to elastic deformation as shown in Table 4. Finally, the validation of code is conducted for different heat source/sink parameter, Q and large values of Prandtl number, Pr, with Liu [66] and Chen [45] keeping k1 = 1, Ec = 0.2, f w = 0 in Table 5. Very good agreement is found in the comparison with minimal percentage errors. Overall therefore confidence in the present hp-FEM computations is very high. For solving above boundary value problem and to give a better approximation for the solution, the suitable guess value of η∞ (length of the domain) is chosen satisfying all boundary conditions. We take the series of values for | − θ  (0) | and | − φ  (0) | with different values of η∞ (such as = 4, 6, 8) with different numbers of elements, E and orders of polynomial, p chosen so that the numerical results obtained are independent of η∞ (see Table 1). For computational purposes, the region of integration η is considered as 0 to η∞ = 6, where η∞ corresponds to η → ∞ which lies significantly outside the momentum and thermal boundary layers.

123

Int. J. Appl. Comput. Math Table 4 Comparison of −θ  (0) among Nataraja et al. [64], Mushtaq et al. [65], Chen [45] and the present results for the PST case with Ec = 0, Q = 0, k1 = 0, f w = 0, N b = N t = 10−5 and no work due to elastic deformation Pr

Nataraja et al. [64]

Mushtaq et al. [65]

Chen [45] (a)

Present results (b)

Percentage error | (b − a) /a| × 100

1

1.3333

1.3349

1.33333

1.33330

0.0018

5

3.3165

3.2927

3.31684

3.31612

0.0218

10

4.7969

4.7742

4.79687

4.79634

0.0110

15

5.9320

5.9097

5.93201

5.93130

0.0120

100

15.7120

15.6884

15.7120

15.70809

0.0249

400

31.6990

31.6289

31.6705

31.65534

0.0478

Table 5 Comparison of −θ  (0) for a second-grade fluid with k1 = 1, Ec = 0.2, f w = 0. Q −0.1

0.0

Pr

Chen [45] (a)

Present results (b)

Percentage error | (b − a) /a| × 100

1

1.37488

1.37488

1.37471

0.0123

10

4.59962

4.59962

4.59893

0.0150

100

14.6843

14.6843

14.6809

0.0231

500

32.8796

32.8796

32.8590

0.0626

1 10

0.1

Liu [66]

–

–

1.34313

–

4.48696

4.48696

4.48601

0.0211

100

14.3328

14.3328

14.3280

0.0335

500

32.0931

32.0931

32.0798

0.0414

1

1.29111

1.29111

1.29109 4.37016

0.0154

10

4.37115

4.37115

100

13.9715

13.9715

13.9621

0.0673

0.0226

500

31.2848

31.2848

31.2677

0.0546

The entire flow domain contains 4001 grid points. At each node four functions are to be evaluated; hence after assembly of the element equations, we obtain a system of 16004 equations which are non-linear. Therefore, an iterative scheme has been employed in the solution. The system is linearized by incorporating the functions f , h and θ , which are assumed to have some prescribed value. After imposing the boundary conditions, a system of 15,097 equations are produced and these are solved by the Gauss elimination method sustaining throughout the computational process an accuracy of 10−4 . The iterative process is terminated when the following condition is satisfied:  −4 |i,m j − i,m−1 (46) j | ≤ 10 i, j

where  denotes either f, h, θ or φ, and m denotes the iterative step. Gaussian quadrature is implemented for solving the integrations. Excellent convergence has been achieved for all the results.

123

Int. J. Appl. Comput. Math

Results and Discussion To provide a physical insight into the present bio-nanopolymer manufacturing flow problem, comprehensive numerical computations are conducted for various values of the parameters that describe the flow characteristics and the results are illustrated graphically. Selected computations are presented in Figs. 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 and 13. In all cases, default values of the governing parameters are: k1 = 0.5, Le = 10, N b = N t = 0.3, Pr = 10, Ec = 0.1, Q = 0.5, f w = 0 unless otherwise stated. These physically correspond to strong viscoelasticity, strong Brownian motion and thermophoresis, weak viscous heating, heat source presence and a solid sheet case (no transpiration at the wall).

Fig. 2 Profiles of stream, velocity, temperature and nanoparticle concentration function for N t = 0.3, N b = 0.3, Pr = 10.0, Le = 10.0, k1 = 0.5, Ec = 0.1, f w = 0, Q = 0.5

Fig. 3 Effect of viscoelastic parameter (k1 ) on stream function distribution with Nt = Nb = 0.3, Pr = 10.0, Le = 10.0, Ec = 0.1, f w = 0, Q = 0.5

123

Int. J. Appl. Comput. Math

Fig. 4 Effect of viscoelastic parameter (k1 ) on velocity distribution with Nt = Nb = 0.3, Pr = 10.0, Le = 10.0, Ec = 0.1, f w = 0, Q = 0.5

Fig. 5 Effect of suction/injection parameter ( f w ) on velocity distribution with Nt = Nb = 0.3, Pr = 10.0, Le = 10.0, Ec = 0.1, k1 = 1.0, Q = 0.05

Figure 2 shows the profiles of stream function ( f ), velocity ( f  ), temperature (θ ) and nanoparticle concentration (φ) for default values of the thermophysical parameters. Smooth profiles are achieved in all cases demonstrating excellent convergence of the finite element computations. Stream function clearly ascends with distance into the boundary layer, whereas velocity, temperature and nanoparticle concentration all descend i.e. these functions are maximized at the wall. Figures 3 and 4 represents the stream function and velocity profiles for different values of viscoelastic parameter, k1 ranging 0–2. It is noted that k1 = 0 is for viscous fluid, k1 > 0 stands for second-grade nanofluid. Significant Brownian motion and thermophoresis are present. The stream function is found to be strongly enhanced with increasing viscoelasticity parameter. All profiles ascend exponentially from zero at the wall to a maximum in the freestream. Fluid velocity however decreases exponentially from unity at the wall to zero at the free stream. Increasing viscoelasticity also elevates the fluid velocity i.e. enhances momentum boundary layer thickness. The viscoelastic nature of the bio-nanofluid therefore benefits the flow and induces acceleration in the boundary layer regime. This trend has been confirmed

123

Int. J. Appl. Comput. Math

Fig. 6 Effect of suction/injection parameter ( f w ) on temperature distribution with Nt = Nb = 0.3, Pr = 10.0, Le = 10.0, Ec = 0.1, k1 = 0.5, Q = 0.05

Fig. 7 Effect of internal heat source/sink parameter (Q) on temperature distribution with Nt = Nb = 0.3, Pr = 10.0, Le = 10.0, Ec = 0.1, k1 = 1.0, f w = 0

in other studies using other viscoelastic non-Newtonian formulations for nanofluids, for example Krishnamurthy et al. [47]. Similar observations have been documented with Jefferys viscoelastic fluid model by Hussain et al. [48] and also the Oldroyd-B model by Khan et al. [49]. Of course these models have a different formulation to the one studied in the current paper; however they do demonstrate similar rheological effects, confirming that the computations elaborated in the present work are in general consistent with other studies. The effects of suction/injection on the velocity and temperature distribution are illustrated in Figs. 5 and 6 respectively, for a second grade nanofluid. As compared to an impermeable sheet ( f w = 0), it is clear that suction ( f w > 0) has the effect to reduce the boundary layer thickness and thus the velocity, whereas injection ( f w < 0) tends to thicken the boundary layer and the velocity increases accordingly. Thus suction acts as a powerful control mechanism for the boundary layer flow i.e. decelerates the flow. Temperature (Fig 6) is also observed to be significantly decreased with increasing suction whereas the converse effect is sustained for increasing injection. Blowing of nanofluid into the boundary layer regime (injection) therefore heats the boundary layer significantly in addition to accelerating the

123

Int. J. Appl. Comput. Math

Fig. 8 Effect of Prandtl number (Pr) on temperature distribution for both (i) Nt = Nb = 10−5 , and (ii) Nt = Nb = 0.3, Le = 10.0, Ec = 0.1, k1 = 1.0, Q = 0.5, f w = 0

Fig. 9 Effect of Eckert number (Ec) on temperature distribution with/without deformation effect keeping Nb = Nt = 0.3, k1 = 1.0, Pr = Le = 10, f w = 0.0, Q = 0.5

flow. Thermal boundary layer thickness is therefore accentuated with an increase of injection with the reverse effect induced with suction (see Fig. 6) In Fig. 7, the effects of temperature dependent heat source/sink (Q) on temperature distribution are shown. The term q (T∞ − T ) signifies the amount of heat generated / absorbed per unit volume, q is a constant, which may take on either positive or negative values. When the wall temperature Tw exceeds the free stream temperature, T∞ , a heat source corresponds to Q > 0 and a heat sink to Q < 0 whereas when Tw < T∞ , the opposite relationship is true. The presence of heat source in the boundary layer generates energy which assists thermal convection and boosts temperatures. This increase in temperature simultaneously accelerates the flow field due to the buoyancy effect. On the other hand, the presence of a

123

Int. J. Appl. Comput. Math

Fig. 10 Variation of heat transfer rate as function of Ec for various values of Nb and Nt keeping Pr = 10.0, Le = 10.0, f w = 0.0

Fig. 11 Effect of Lewis number (Le) on temperature distribution with Nb = Nt = 0.1, Ec = 0.1, k1 = 1.0, f w = 0, Q = 0.5

heat sink in the boundary layer absorbs energy which causes the temperature of the fluid to decrease. Thermal boundary layer thickness of the viscoelastic biopolymer nanofluid sheet will be increased with a heat source and depleted with a heat sink. The effects of Brownian motion parameter, Nb and thermophoresis parameter, Nt, on temperature are shown in Fig. 8. As expected, the boundary layer profiles for the temperature are of the same form as in the case of regular viscoelastic fluids. The temperature in the boundary layer increases with the increase in the Brownian motion parameter (Nb) and thermophoresis parameter (Nt). The Brownian motion of nanoparticles can enhance thermal conduction via several methods including for example, direct heat transfer owing to nanoparticles or by virtue of micro-convection of fluid surrounding individual nanoparticles. For larger diameter particles, Brownian motion will be weaker and the parameter, Nb will have lower values. For

123

Int. J. Appl. Comput. Math

Fig. 12 Effect of Lewis number (Le) on nanoparticle concentration with Nb = Nt = 0.1, Ec = 0.1, k1 = 1.0, f w = 0, Q = 0.5

Fig. 13 Variation of mass transfer rate as function of Q for various suction/injection parameter keeping Le = 10.0, Nb = Nt = 0.3, Ec = 0.1

smaller diameter particles Brownian motion will be greater and Nb will have larger values. In accordance with this, we observe that temperatures are enhanced with higher Nb values whereas they are reduced with lower Nb values. Brownian motion therefore contributes significantly to thermal enhancement in the boundary layer regime (Fig. 8). Similarly increasing thermophoresis (Nt) which is due to temperature gradient and associated with particle deposition, also leads to an increase in the temperature profile, as witnessed in Fig. 8. Furthermore Fig. 8 also exhibits the reduction in temperatures caused by an increase in Prandtl number. The larger values of Prandtl number (Pr) imply a much lower thermal conductivity of the viscoelastic bio-nanofluid which serves to depress thermal diffusion and cools the boundary layer regime.

123

Int. J. Appl. Comput. Math

Figure 9 illustrates the response of temperature profiles to a variation in Eckert number with/without work done due to deformation keeping Nb = Nt=0.5, k1 = 0.5, Pr = Le = 10, f w = 0.1, Q = 1.0. Viscous heating enhances temperatures and thickens the thermal boundary layer. However the increase is markedly more pronounced for the case of work done due to deformation, rather than in absence of work done due to deformation, for high value of Eckert number. Figure 10 presents the variation in dimensionless heat transfer rates with Eckert number, and furthermore includes the influence of Nb and Nt parameters on the dimensionless heat transfer rates. Viscous dissipation (as characterized by the Eckert number) and work done by deformation strongly decrease the heat transfer, since greater thermal energy is dissipated in the boundary layer regime and this results in a depletion of heat transferred to the wall. Moreover, heat transfer rate is also decreased with the increase of Brownian motion and thermophoresis, since as established earlier both Brownian motion and thermophoresis enhance boundary layer temperatures leading to a reduction in transport of heat to the wall. These trends concur with the earlier computations of Khan and Pop [26]. It is evident overall from Fig. 10 that the dimensionless heat transfer rate is a decreasing function of Nb, Nt and Ec. Figures 11 and 12 depict the variation of temperature and nanoparticle concentration for various Lewis numbers (Le). Lewis number defines the ratio of thermal diffusivity to mass diffusivity. It is used to characterize fluid flows where there is simultaneous heat and mass transfer by convection. Effectively, it is also the ratio of Schmidt number and the Prandtl number. Temperature and thermal boundary layer thickness are slightly decreased with an increase in Lewis number (Fig. 11). Nanoparticle concentration function, φ (η), is however found to be very significantly reduced with increasing Lewis number (Fig. 12). This is attributable to the decrease in mass (species) diffusivity associated with an increase in Lewis number. Species diffusion rate is therefore depressed as Lewis number increases which manifests in a strong fall in concentrations. 1/2 Figure 13 depicts the distributions of the mass transfer function (Sh x Rex ) with heat source/sink parameter (Q) for different values of suction/injection parameter. The mass transfer increases with increase of heat source (Q > 0) whereas it is decreased with increasing heat sink parameter (Q < 0). An increase in injection parameter ( f w < 0) strongly suppresses the mass transfer at the wall whereas increasing suction is found to enhance it. The presence of a heat source and wall suction therefore have significant beneficial effects on transport phenomena in stretching sheet nanofluid processing, whereas a strong heat sink and blowing (injection) tend to inhibit transport.

Conclusions In the present paper, a mathematical model is developed for viscoelastic bio-nano-polymer extrusion from a stretching sheet with Brownian motion and thermophoresis effects incorporated. The governing partial differential equations for mass, momentum, energy and species conservation are rendered into a system of coupled, nonlinear, ordinary differential equations by using a similarity transformation. The higher order finite element method (hp-FEM) has been implemented to solve the resulting two-point nonlinear boundary value problem more efficiently. Excellent correlation with previous published results has been achieved. The computations have shown that: 1. An increase in the polymer fluid viscoelasticity (k1 ) accelerates the flow.

123

Int. J. Appl. Comput. Math

2. Increasing Brownian motion parameter (Nb) and thermophoresis parameter (Nt) enhance temperature in the boundary layer region whereas they reduce the heat transfer rates (local Nusselt number function). 3. The kinetic energy dissipation (represented by the Eckert number, Ec) due to viscous heating and deformation work has the effect to thicken the thermal boundary layer and strongly elevates temperatures in the viscoelastic nano-bio-polymer. 4. Increasing the Lewis number (Le) decreases temperature weakly whereas it strongly reduces nanoparticle concentrations. 5. An increase in Prandtl number (Pr) significantly decreases temperatures. 6. The presence of internal heat generation (Q > 0) enhances temperatures and therefore reduces the heat transfer rate (local Nusselt number function), with the opposite trend sustained for the case of heat absorption (Q < 0) for nanofluid. 7. Increasing suction ( f w > 0) strongly decelerates the nanofluid boundary layer flow, decreases nanofluid temperatures and enhances mass transfer rates (local Sherwood number function), whereas increasing injection ( f w < 0) accelerates the flow, enhances temperatures and depresses wall mass transfer rates. The present hp-FEM shows excellent accuracy and stability and will be employed in further simulating flows of interest in bio-nano-polymer manufacturing processes involving other viscoelastic models e.g. Maxwell fluids [67] and also nano-particle geometry effects [68]. Acknowledgments Dr. O. Anwar Bég is grateful to the late Professor Howard Brenner (1929-2014) of Chemical Engineering, MIT, USA, for some excellent discussions regarding viscoelastic characteristics of biopolymers.

References 1. Nie, Z., Kumacheva, E.: Patterning surfaces with functional polymers. Nat. Mater. 7, 277–290 (2008) 2. Pal, K., Banthia, A.K., Majumdar, D.K.: Polymeric hydrogels: characterization and biomedical applications. Des. Monomers Polym. 12, 197–220 (2009) 3. Becker, T.A., Kipke, D.R.: Flow properties of liquid calcium alginate polymer injected through medical microcatheters for endovascular embolization. J. Biomed. Mater. Res. 61, 533–540 (2002) 4. Richard, A., Margantis, A.: Production and mass transfer characteristics of non-Newtonian biopolymers for biomedical applications. Crit. Rev. Biotechnol. 22, 355–374 (2002) 5. Choi, S.: Enhancing thermal conductivity of fluids with nanoparticles. In: Siginer, D.A., Wang, H.P. (eds.) Developments and applications of non-Newtonian flows, vol. 66, pp. 99–105. ASME Fluids Engineering Division, Boston (1995) 6. Anwar Bég, O., Tripathi, D.: Mathematica simulation of peristaltic pumping with double-diffusive convection in nanofluids: a bio–nano-engineering model. Proc. Inst. Mech. Eng. Part N J. Nanoeng. Nanosyst. 225, 99–114 (2012) 7. Anwar Bég, O., Rashidi, M.M., Akbari, M., Hosseini, A.: Comparative numerical study of single-phase and two-phase models for bio-nanofluid transport phenomena. J. Mech. Med. Biol. 14, 1450011.1– 1450011.31 (2014) 8. Pak, B.C., Cho, Y.: Hydrodynamics and heat transfer study of dispersed fluids with submicron metallic oxide particles. Exp. Heat Transf. 11, 151–170 (1998) 9. Karthikeyan, N.R., Philip, J., Raj, B.: Effect of clustering on the thermal conductivity of nanofluids. Mater. Chem. Phys. 109, 50–55 (2008) 10. Prakash, M., Giannelis, E.P.: Mechanism of heat transport in nanofluids. J. Comput. Aided Mater. Des. 14, 109–117 (2007) 11. Wang, X., Xu, X., Choi, S.U.S.: Thermal conductivity of nanoparticle fluid mixture. AIAA J. Thermophys. Heat Trans. 13(4), 474–480 (1999) 12. Jang, S.P., Choi, S.U.S.: Role of Brownian motion in the enhanced thermal conductivity of nanofluids. App. Phys. Lett. 84, 4316–4318 (2004)

123

Int. J. Appl. Comput. Math 13. Chon, C.H., Kihm, K.D., Lee, S.P., Choi, S.U.S.: Empirical correlation finding the role of temperature and particle size for nanofluid (Al2 O3 ) thermal conductivity enhancement. App. Phys. Lett. 87, 153107 (2005) 14. Keblinski, P., Eastman, J.A., Cahill, D.G.: Nanofluids for thermal transport. Mat. Today 8, 36–44 (2005) 15. Leong, K.C., Yang, C., Murshed, S.M.S.: A model for the thermal conductivity of nanofluids: the effect of interfacial layer. J. Nanoparticle Res. 8, 245–254 (2006) 16. Buongiorno, J.: Convective transport in nanofluids. ASME J. Heat Transf. 128, 240–250 (2006) 17. Eastman, J.A., Phillpot, S.R., Choi, S.U.S., Keblinski, P.: Thermal transport in nanofluids. Ann. Rev. Mater. Res. 34, 219–146 (2004) 18. Wang, X.-Q., Majumdar, A.S.: Heat transfer characteristics of nanofluids: a review. Int. J. Thermal Sci. 46, 1–19 (2007) 19. Trisaksri, V., Wongwises, S.: Critical review of heat transfer characteristics of nanofluids. Renew. Sustain. Energy Rev. 11, 512–523 (2007) 20. Wen, D., Lin, G., Vafaei, S., Zhang, K.: Review of nanofluids for heat transfer applications. Particuology 7, 141–150 (2009) 21. Saidur, R., Leong, K.Y., Mohammad, H.A.: A review on applications and challenges of nanofluids. Renew. Sustain. Energy Rev. 15, 1646–1668 (2011) 22. Kuznetsov, A.V., Nield, D.A.: Natural convection boundary layer flow of a nanofluids past a vertical plate. Int. J. Therm. Sci. 49, 243–247 (2010) 23. Rana, P., Bhargava, R., Anwar Bég, O.: Numerical solution for mixed convection boundary layer flow of a nanofluid along an inclined plate embedded in a porous medium. Comput. Math. Appl. 64, 2816–2832 (2012) 24. Rana, P., Bhargava, R.: Flow and heat transfer analysis of a nanofluid along a vertical flat plate with non-uniform heating using FEM: effect of nanoparticle diameter. Int. J. Appl. Phys. Math. 1, 171–176 (2011) 25. Prasad, V.R., Gaffar, S.A., Anwar Bég, O.: Heat and mass transfer of a nanofluid from a horizontal cylinder to a micropolar fluid. AIAA J. Thermophys. Heat Transf. 29, 127–139 (2015) 26. Khan, W.A., Pop, I.: Boundary-layer flow of a nanofluid past a stretching sheet. Int. J. Heat Mass Transf. 53, 2477–2483 (2010) 27. Rana, P., Bhargava, R.: Flow and heat transfer of a nanofluid over a nonlinearly stretching sheet: a numerical study. Commun. Nonlinear Sci. Num. Simul. 17(1), 212–226 (2012) 28. Uddin, M.J., Anwar Bég, O., Amin, N.S.: Hydromagnetic transport phenomena from a stretching or shrinking nonlinear nanomaterial sheet with Navier slip and convective heating: a model for bio-nanomaterials processing. J. Magnetism Magn. Mater. 368, 252–261 (2014) 29. Uddin, M.J., Bég, O.A., Ismail, A.I.M: Mathematical modelling of radiative hydromagnetic thermosolutal nanofluid convection slip flow in saturated porous media. Math. Prob. Eng. Article ID 179172, 11 pages. doi:10.1155/2014/179172 (2014) 30. Nadeem, S., Lee, C.: Boundary layer flow of nanofluid over an exponentially stretching surface. Nanoscale Res. Lett. 7, 94 (2012) 31. Kandasamy, R., Loganathan, P., PuviArasu, P.: Scaling group transformation for MHD boundary-layer flow of a nanofluid past a vertical stretching surface in the presence of suction/injection. Nuclear Eng. Des. 241(6), 2053–2059 (2011) 32. Bachok, N., Ishak, A., Pop, I.: Unsteady boundary-layer flow and heat transfer of a nanofluid over a permeable stretching/shrinking sheet. Int. J. Heat Mass Transf. 55(7–8), 2102–2109 (2012) 33. Rana, P., Bhargava, R., Anwar Bég, O.: Unsteady MHD transport phenomena on a stretching sheet in a rotating nanofluid. Proc. Int. Mech. E Part N J. Nanoeng. Nanosyst. 227, 277–299 (2013) 34. Chen, H.S., Ding, T.L., He, Y.R., Tan, C.Q.: Rheological behaviour of ethylene glycol based titanium nanofluids. Chem. Phys. Lett. 444, 333–337 (2007) 35. Chen, H.S., Ding, T.L., Lapkin, A.: Rheological behaviour of nanofluids containing tube/rod like nanoparticles. Power Technol. 194, 132–141 (2009) 36. Gallego, M.J.P., Lugo, L., Legido, J.L., Piñeiro, M.M.: Rheological non-Newtonian behaviour of ethylene glycol-based Fe2 O3 -nanofluids. Nanoscale Res. Lett. 6(1), 560 (2011) 37. Khan, W.A., Gorla, R.S.R.: Heat and mass transfer in non-Newtonian nanofluids over a non-isothermal stretching wall. Proc. Inst. Mech. E Part N J. Nanoeng. Nanosyst. 225, 155–163 (2011) 38. Rajagopal, K.R., Na, T.Y., Gupta, A.S.: Flow of a viscoelastic fluid over a stretching sheet. Rheol. Acta 23, 213–215 (1984) 39. Rajagopal, K.R., Na, T.Y., Gupta, A.S.: A non-similar boundary layer on a stretching sheet in a nonNewtonian fluid with uniform free stream. J. Math. Phys. Sci. 21(2), 189–200 (1987) 40. Dandapat, B.S., Gupta, A.S.: Flow and heat transfer in a viscoelastic fluid over a stretching sheet. Int. J. Nonlinear Mech. 24(3), 215–219 (1989)

123

Int. J. Appl. Comput. Math 41. Rao, B.N.: Technical note: flow of a fluid of second grade over a stretching sheet. Int. J. Nonlinear Mech. 31(4), 547–550 (1996) 42. Khan, S.K.: Heat transfer in a viscoelastic fluid flow over a stretching surface with heat source/sink, suction/blowing and radiation. Int. J. Heat Mass Transf. 49(3–4), 628–639 (2006) 43. Anwar Bég, O., Bég, T.A., Takhar, H.S., Raptis, A.: Mathematical and numerical modeling of nonNewtonian thermo-hydrodynamic flow in non-Darcy porous media. Int. J Fluid Mech. Res. 31, 1–12 (2004) 44. Anwar Bég, O., Takhar, H.S., Bharagava, R., Rawat, S., Prasad, V.R.: Numerical study of heat transfer of a third grade viscoelastic fluid in non-Darcian porous media with thermophysical effects. Phys. Scr. 77, 1–11 (2008) 45. Chen, C.-H.: On the analytic solution of MHD flow and heat transfer for two types of viscoelastic fluid over a stretching sheet with energy dissipation, internal heat source and thermal radiation. Int. J. Heat Mass Transf. 53, 4264–4273 (2010) 46. Anwar Bég, O., Sharma, S., Bhargava, R., Bég, T.A.: Finite element modelling of transpiring third-grade viscoelastic biotechnological fluid flow in a Darcian permeable half-space. Int. J. Appl. Math. Mech. 7, 38–52 (2011) 47. Krishnamurthy, M.R., Prasannakumara, B.C., Gireesha, B.J., Gorla, R.S.R.: Effect of chemical reaction on MHD boundary layer flow and melting heat transfer of Williamson nanofluid in porous medium. Eng. Sci. Technol. Int. J. 19, 53–61 (2016) 48. Hussain, T., Shehzad, S.A., Hayat, T., Alsaedi, A., Al-Solamy, F.: Radiative hydromagnetic flow of Jeffrey nanofluid by an exponentially stretching sheet. PLoS ONE 9(8), e103719 (2014). doi:10.1371/journal. pone.0103719 49. Khan, W.A., Khan, M., Malik, R.: Three-dimensional flow of an Oldroyd-B nanofluid towards stretching surface with heat generation/absorption. PLoS ONE 9(8), e105107 (2014). doi:10.1371/journal.pone. 0105107 50. Sher Akbar, N., Ebaid, A., Khan, Z.H.: Numerical analysis of magnetic field on Eyring–Powell fluid flow towards a stretching sheet. J. Magn. Magn. Mater. 382, 355–358 (2015) 51. Mehmood, R., Nadeem, S., Akbar, N.S.: Oblique stagnation flow of Jeffery fluid over a stretching convective surface: optimal solution. Int. J. Numer. Methods Heat Fluid Flow 25(3), 454–471 (2015) 52. Nadeem, S., Akbar, N.S., Khan, Z.H.: Buoyancy and radiation effect on stagnation point flow of micropolar nanofluid along a vertically convective stretching surface. IEEE Trans. Nanotechnol. 14, 42–50 (2015) 53. Babuska, I., Guo, B.Q.: The h, p and h–p version of the finite element method: basic theory and applications. Advs. Eng. Softw. 15, 159–174 (1992) 54. Khomami, B., Talwar, K.K., Ganpule, H.K.: A comparative study of higher-and lower-order finite element techniques for computation of viscoelastic flows. J. Rheol. 38, 255–289 (1994) 55. Coleman, B.D., Noll, W.: An approximation theorem for functionals with applications in continuous mechanics. Arch. Ration. Mech. Anal. 6, 355–370 (1960) 56. Dunn, J.E., Fosdick, R.L.: Thermodynamics, stability and boundedness of fluids of complexity 2 and fluids of second grade. Arch. Ration. Mech. Anal. 56, 191–252 (1974) 57. Fosdick, R.L., Rajagopal, K.R.: Anomalous features in the model of second order fluids. Arch. Ration. Mech. Anal. 70, 145–152 (1979) 58. Dunn, J.E., Rajagopal, K.R.: Fluids of differential type: critical review and thermodynamic analysis. Int. J. Eng. Sci. 33, 689–729 (1995) 59. Troy, W.C., Overman II, E.A., Ermentrout, G.B., Keener, J.P.: Uniqueness of a second order fluid past a stretching sheet. Quart. Appl. Math. 45, 755–793 (1987) 60. McLeod, B., Rajagopal, K.R.: On uniqueness of flow of a Navier–Stokes fluid due to stretching boundary. Arch. Mech. Anal. 98, 393–985 (1987) 61. Chang, W.D.: The non-uniqueness of the flow of viscoelastic fluid over a stretching sheet. Quart. Appl. Math. 47, 365–366 (1989) 62. Wang, C.Y.: Free convection on a vertical stretching surface. J. Appl. Math. Mech. (ZAMM) 69, 418–420 (1989) 63. Gorla, R.S.R., Sidawi, I.: Free convection on a vertical stretching surface with suction and blowing. Appl. Sci. Res. 52, 247–257 (1994) 64. Nataraja, H.R., Sarma, M.S., Rao, B.N.: Flow of a second-order fluid over a stretching surface having power-law temperature. Acta Mech. 128, 259–262 (1998) 65. Mushtaq, M., Asghar, S., Hossain, M.A.: Mixed convection flow of second grade fluid along a vertical stretching flat surface with variable surface temperature. Heat Mass Transf. 43, 1049–1061 (2007) 66. Liu, I.-C.: Flow and heat transfer of an electrically-conducting fluid of second grade over a stretching sheet subject to a transverse magnetic field. Int. J. Heat Mass Transf. 47, 4427–4437 (2004)

123

Int. J. Appl. Comput. Math 67. Norouzi, M., Davoodi, M., Anwar Bég, O.: An analytical solution for convective heat transfer of viscoelastic flows in rotating curved pipes. Int. J. Therm. Sci. 90, 90–111 (2015) 68. Rana, P., Anwar Bég, O.: Mixed convection flow along an inclined permeable plate: effect of magnetic field, nanolayer conductivity and nanoparticle diameter. Appl. Nanosci. 5(5), 569–581 (2015)

123