NUMERICAL STUDY OF NON-NEWTONIAN POLYMERIC BOUNDARY LAYER FLOW AND HEAT TRANSFER FROM A PERMEABLE HORIZONTAL ISOTHERMAL CYLINDER

In this article, we investigate the nonlinear steady state boundary layer flow and heat transfer of an incompressible Jeffery non-Newtonian fluid from a permeable horizontal isothermal cylinder. The transformed conservation equations are solved numerically subject to physically appropriate boundary conditions using a versatile, implicit, finite-difference technique. The numerical code is validated with previous studies. The influence of a number of emerging non-dimensional parameters, namely with Deborah number (De), surface suction parameter (S), Prandtl number (Pr), ratio of relaxation to retardation times (λ) and dimensionless tangential coordinate (ξ) on velocity and temperature evolution in the boundary layer regime are examined in detail. Furthermore, the effects of these parameters on surface heat transfer rate and local skin friction are also investigated. It is found that the velocity is reduced with increasing Deborah number whereas temperature is enhanced. Increasing λ enhances the velocity but reduces the temperature. The heat transfer rates is found to be depressed with increasing Deborah number, De, and enhanced with increasing λ. Local skin friction is found to be decreased with a rise in Deborah number whereas it is elevated with increasing values of relaxation to retardation time ratio (λ). Increasing suction decelerates the flow and cools the boundary layer i.e. reduces temperatures. With increasing tangential coordinate, the flow is also decelerated whereas the temperatures are enhanced. The simulation is relevant to polymer coating thermal processing. Polymeric enrobing flows are important in industrial manufacturing technology and process systems. Such flows are non-Newtonian. Motivated by such applications, we did the present problem.


INTRODUCTION
Polymeric flows are generally non-Newtonian in nature. In coating applications, fluid mechanics and heat transfer play a key role in determining the constitution of manufactured polymers Middleman (1997). In dip coating processes for example the surface to be coated is first immersed in a polymer and then steadily withdrawn Roy (1971). Since polymers have high viscosity levels, industrial chemical/manufacturing flow processes exploiting such fluids are generally laminar in nature Baaijens et.al (1997). The rheological nature of polymers also necessitates more sophisticated mathematical models for describing shear stress-strain relationships Otsuki (1999). Important characteristics include relaxation, retardation, viscoelasticity and elongational viscosity. A number of theoretical and computational studies have been communicated on transport phenomena from cylindrical bodies, which frequently arise in polymer processing systems. These Newtonian studies were focused more on heat transfer aspects and include Eswara and Nath Eswara and Nath (1992), Rotte and Beek (1962), and the pioneering analysis of Sakiadis (1961). Further more recent studies examining multi-physical and chemical transport from cylindrical bodies include Zueco et al. (2009Zueco et al. ( , 2011). An early investigation of rheological boundary layer heat transfer from a horizontal cylinder was presented by Chen and Leonard (1972), who considered the power-law model and demonstrated that the transverse curvature has a strong effect on skin friction at moderate and large distances from the leading edge of the boundary layer. Lin and Chen (1979) also studied axisymmetric laminar boundary-layer convection flow of a power-law non-Newtonian fluid over both a circular cylinder and a spherical body using the Merk-Chao series solution method. Pop et al. (1990) simulated numerically the steady laminar forced convection boundary layer of power-law non-Newtonian fluids on a continuously moving cylinder with the surface maintained at a uniform temperature or uniform heat flux. Further non-Newtonian models employed in analyzing convection flows from cylinders include micropolar liquids chang (2006), viscoelastic materials Anwar et.al (2008) and Kasim et.al (2011), micro polar nanofluids Rehman and Nadeem (2011) and Casson fluids Prasad et.al (2012). One subclass of non-Newtonian fluids known as the Jeffery fluid Saasen and Hassager (1991) is particularly useful owing to its simplicity. This fluid model is capable of describing the characteristics of relaxation and retardation times, which arise in complex polymeric flows. Furthermore, the Jeffrey type model utilizes time derivatives rather than convected derivatives, which make it more amenable for numerical simulations. Recently the Jeffery model has received considerable attention.

2
The objective of the present paper is to investigate the laminar boundary layer flow and heat transfer of a Jeffery non-Newtonian fluid from a horizontal cylinder. The non-dimensional transport equations with associated dimensionless boundary conditions constitute a highly nonlinear, coupled two-point boundary value problem. Keller's implicit finite difference "box" scheme is implemented to solve the problem. The effects of the emerging thermophysical parameters, namely Prandtl number, Deborah number, ratio of relaxation to retardation times and transpiration (wall suction or injection) on the momentum and heat transfer characteristics are studied. The present problem has to the authors' knowledge not appeared thus far in the scientific literature and is relevant to polymeric manufacturing processes.

MATHEMATICAL MODEL
Laminar incompressible boundary layer flow and heat transfer under thermal buoyancy force, from a horizontal permeable cylinder to a Jeffery rheological fluid is examined, as illustrated in Figure 1. Both the cylinder and the Jeffery fluid are maintained initially at the same temperature. Instantaneously they are raised to a temperature w where the latter (ambient) temperature of the fluid is sustained constant. The x-coordinate (tangential) is orientated along the circumference of the horizontal cylinder from the lowest point and the y-coordinate (radial) is directed perpendicular to the surface, with a denoting the radius of the horizontal cylinder. Ф=x/a represents the angle of the yaxis with respect to the vertical (0≤Ф≤π). The gravitational acceleration g, acts vertically downwards. The Boussineq approximation holds i.e. density variation is only experienced in the buoyancy term in the momentum equation. The Cauchy stress tensor, S, of a Jeffrey non-Newtonian fluid [24] takes the form: where a dot above a quantity denotes the material time derivative, p is pressure, I is the identity tensor,  is dynamic viscosity, λ is the ratio of relaxation to retardation times, λ1 is the retardation time and  is the shear rate. The Jeffery model provides an elegant formulation for simulating retardation and relaxation effects arising in polymer flows. Introducing the boundary layer approximations, and incorporating the stress tensor for a Jeffery fluid in the momentum equation (in differential form) the conservations equations take the form: x y x y y x y y y The Jeffery fluid model therefore introduces a number of mixed derivatives into the momentum boundary layer equation (3) and in particular two third order derivatives, making the system an order higher than the classical Navier-Stokes (Newtonian) viscous flow model. The non-Newtonian effects feature in the shear terms only of eqn.
(3) and not the convective (acceleration) terms. The final term on the right hand side of eqn.
(3) represents the thermal buoyancy effect (free convection) and couples eqn.
(3) to the energy equation (4). Viscous dissipation effects are neglected in the model. In eqns.
(2)-(4) u and v designate velocity components in the x -and y -directions respectively,     is the kinematic viscosity of the Jeffery fluid,  is the coefficient of thermal expansion, is the thermal diffusivity, T is the temperature, k is the thermal conductivity of the Jeffery fluid, ρ is the density of the Jeffery fluid, cp is the specific heat at constant pressure, and all other parameters have been defined earlier. The appropriate boundary conditions are imposed at the cylinder surface and in the free stream (edge of the boundary layer) and take the form: At , and therefore, the mass conservation eqn.
(2) is automatically satisfied. Furthermore, the following dimensionless variables are introduced into eqns.
(2)-(4) and (5): , Pr Where  is the non-dimensional tangential coordinate, is the nondimensional radial coordinate, f is dimensionless stream function, Pr is the Prandtl number, S is the wall transpiration parameter (S> 0 for suction and S< 0 for injection, the case of an impermeable cylinder is retrieved for S = 0), Gr is the Grashof number and De is the Deborah number characterizing the fluidity of the material (viscoelasticity). The resulting momentum and thermal boundary layer equations take the form: The location,  0, corresponds to the vicinity of the lower stagnation point on the cylinder, an aspect discussed in more detail subsequently. The corresponding non-dimensional boundary conditions for the collectively sixth order, multi-degree partial differential equation system defined by eqns. (7), (8) assume the form: Here primes denote the differentiation with respect to .The skinfriction coefficient (shear stress at the cylinder surface) and Nusselt number (heat transfer rate at the cylinder surface) can be defined using the transformations described above with the following expressions:

COMPUTATIONAL SOLUTION WITH KELLER BOX IMPLICT METHOD
In Purely analytical solutions for the boundary value problem defined by eqns. (7), (8) and conditions (10) are extremely difficult, if not intractable. A computational solution is therefore developed using a versatile and stable finite difference algorithm based on Keller's box method Keller (1970). This implicit method has been implemented extensively in non-Newtonian fluid mechanics simulations for a variety of different rheological models. Hossain et al. (1970) used a power-law model and Keller's scheme to study free convection boundary layers from a slotted vertical plate. Javed et al. (2013)   . The present code has received extensive validation in previous studies, as described in Subbarao et.al (2016) and therefore confidence is high in the accuracy of computations. The fundamental steps of the Keller Box Scheme are as follows: 1) Reduction of the N th order partial differential equation system to N first order equations 2) Finite Difference Discretization 3) Quasilinearization of Non-Linear Keller Algebraic Equations 4) Block-tridiagonal Elimination of Linear Keller Algebraic Equations.
Step 1: Reduction of the N th order partial differential equation system to N first order equations New variables are introduced to Eqns. (7) - (8) and (9), to render the boundary value problem as a multiple system of first order equations. A set of six simultaneous first order differential equations are therefore generated by introducing the new variables u, v, q and t: where primes denote differentiation with respect to  . In terms of the dependent variables, the boundary conditions become:

Fig.2 Keller Box element and boundary layer mesh
Step 2: Finite Difference Discretization A two dimensional computational grid is imposed on the -η plane as sketched in Fig. 2. The stepping process is defined by:   De vv Where the following notation applies   The boundary conditions are Step

3: Quasilinearization of Non-Linear Keller Algebraic Equations
Assuming f u v q s t , j = 0, 1, 2 …, J. This nonlinear system of algebraic equations is linearized by means of Newton's method as explained in Subba Rao et.al. (2015.

Step 4: Block-tridiagonal Elimination of Linear Keller Algebraic Equations
The linear system (15a) -(15f) can now be solved by the blockelimination method owing to its block-tridiagonal structure. The bocktridiagonal structure generated consists of block matrices. The complete linearized system is formulated as a block matrix system, where each element in the coefficient matrix is a matrix itself, and this system is solved using the efficient Keller-box method. The numerical results are strongly influenced by the number of mesh points in both directions. After some trials in the η-direction (radial coordinate), a larger number of mesh points are selected whereas in the ξ-direction (tangential coordinate) significantly less mesh points are necessary. max  has been set at 10 and this constitutes an adequately large value at which the prescribed boundary conditions are satisfied. max  is set at 3.0 for the simulations. Mesh independence has been comfortably attained in the simulations. The numerical algorithm is executed in MATLAB.

RESULTS AND DISCUSSION
Comprehensive solutions have been obtained and are presented in Tables 1 and 2 and Figs. 3-7. Table 1 presents the influence of increasing Prandtl number on skin friction and heat transfer rate, along with a variation in the parameter (λ), ratio of relaxation and retardation times and transverse coordinate (). With increasing Prandtl number, the skin friction is generally decreased, whereas heat transfer rate is markedly enhanced. Heat transfer rate is maximized at the lower stagnation point ( =0), for any value of Prandtl number or rheological parameter (). With an increase in , both skin friction and heat transfer rates are increased. This implies that as the relaxation time is reduced (and the retardation time increased) the polymer flows faster and transfers heat more efficiently from the cylinder surface. This appears consistent with other studies Hayat et.al. (2012). in table 1, S is positive corresponding to wall suction.  Table 2 presents the influence of increasing Deborah number parameter, Deon skin friction and heat transfer rate, along with a variation in the suction parameter, S and transverse coordinate, . With increasing Deborah number, the skin friction is generally decreased, and heat transfer rate is also decreased. This trend is sustained for all values of transverse coordinate. An increase in suction (S> 0) reduces skin friction at higher values of the transverse coordinate, for any value of Deborah number. Further from the vicinity of the lower stagnation point, therefore the polymer flow is decelerated with suction. However, at lower values of transverse coordinate, suction slightly accelerates the flow for all Deborah numbers. In the vicinity of the lower stagnation point, 0 and the boundary layer equations eqns. (7) to (8) contract to a system of ordinary differential equations:  De It therefore is intimately associated with the shearing characteristics of the polymer flow. For polymers, larger De values imply that the polymer becomes highly oriented in one direction and stretched. Generally, this arises when the polymer takes longer to relax in comparison with the rate at which the flow is deforming it. When such fluids are stretched, there is a delay in their return to the unstressed state. For very large Deborah numbers, the fluid movement is too fast for elastic forces to relax and the material then acts like a purely elastic solid. Large Deborah numbers are therefore not relevant to the present simulations. For small Deborah numbers, the time scale of fluid movement is much greater than the relaxation time of elastic forces in the polymer and the polymer then behaves as a simple viscous fluid, as elaborated by Bég and Makinde (2011). Vrentas et al. (1975) have also indicated that the Deborah number can be utilized in characterizing diffusional transport in amorphous polymer-solvent systems. Further from the cylinder surface, we observe that there is a slight increase in velocity i.e. the flow is accelerated with increasing Deborah number. With greater distance from the solid boundary, the polymer is therefore assisted in flowing even with higher elastic effects. Clearly, the responses in the near-wall region and far-field region are very different. In fig. 3b, an increase in Deborah number is seen to considerably enhance temperatures throughout the boundary layer regime. This has also been observed by Hayat et al. (2012). Although De does not arise in the thermal boundary layer equation (8) ff The thermal boundary layer equation (8) remains unchanged. Effectively with greater relaxation time of the polymer, the thermal boundary layer thickness is reduced. However, with greater relaxation times, the momentum boundary layer thickness is only decreased near the cylinder surface whereas further away it is enhanced since the flow is strongly accelerated there.  6(a) -6(b) depict the velocity and temperature distributions with radial coordinate, for various transverse (stream wise) coordinate values, . Generally, velocity is noticeably lowered with increasing migration from the leading edge i.e. larger  values (figure 6a). The maximum velocity is computed at the lower stagnation point (~0) for low values of radial coordinate (). The transverse coordinate clearly exerts a significant influence on momentum development. A very strong increase in temperature (), as observed in figure 6b, is generated throughout the boundary layer with increasing  values. The temperature field decays monotonically. Temperature is maximized at the cylinder surface (=0) and minimized in the free stream (= 8).
Although the behavior at the upper stagnation point (~) is not computed, the pattern in figure 6b suggests that temperature will continue to progressively grow here compared with previous locations on the cylinder surface (lower values of ). Similar observations have been recorded by a number of researchers including Chen and Leonard (1972) for power-law fluids, Chang (2006) for micro polar fluids and Anwar et al. (2008) for viscoelastic fluids. Figures 7(a) -7(b) present typical profiles for velocity ( f  ) and temperature for various values of the transpiration parameter, S. As in all other graphs, only the case of wall suction is studied (S > 0). It is observed that an increase in the suction parameter significantly decelerates the flow for all values of radial coordinate. The boundary layer thickness is reduced and suction causes the boundary layer to 7 adhere closer to the wall. Similarly increasing wall suction is found to lower temperatures in the boundary layer regime and strongly decreases thermal boundary layer thickness. Although boundary layer separation has not been identified in the present regime, suction has been shown to delay this effect in certain viscoelastic cylinder flow problems.

CONCLUSIONS
A mathematical model has been developed for boundary layer mixed convection flow of a Jeffery non-Newtonian fluid from a horizontal cylinder, with wall transpiration. The transformed conservation equations have been solved with prescribed boundary conditions using the implicit Keller-box finite difference method. The present simulations have shown that: 1. Increasing the viscoelasticity parameter i.e. Deborah number (De), reduces the velocity, skin friction and heat transfer rate whereas it enhances temperature.
2. Increasing the parameter ratio of relaxation and retardation times (λ), increases velocity, skin friction coefficient and heat transfer rate whereas it reduces temperature.
4. Increasing transverse coordinate () generally decelerates the flow near the cylinder surface and reduces momentum boundary layer thickness whereas it enhances temperature. Heat transfer rate is also maximized at the lower stagnation point (= 0). 5. Increasing suction at the cylinder surface (fw>0) decelerates the flow and strongly depresses temperature.
Generally, very stable and accurate solutions are obtained with the present finite difference code and it is envisaged that other non-Newtonian flows will be studied using this methodology in the future including Maxwell upper convected fluids (Vajravelu,2012) and couple stress fluids (Beg et al., 2012).

Cf
skin friction coefficient f non-dimensional steam function g acceleration due to gravity k thermal conductivity of Jeffery fluid cp specific heat capacity Gr Grashof (