Numerical Solutions for Micropolar Transport Phenomena over a Nonlinear Stretching Sheet

Department of Mathematics, Indian Institute of Technology, Roorkee, I ndia Department of Engineering, Manchester Metropolitan University Manchester, M1 5GD, England, UK 3 Fire Safety Engineering Science Program Leeds College of Building, Leeds Metropolitan University North Street, Leeds, LS2 7QT, England obeg@lcb.ac.uk; docoabeg@hotmail.com Department of Civil Engineering, Indian Institute of Technology, Roor kee, India

f dimensionless stream function T temperature g dimensionless microrotation u velocity along the x-axis g e gravitational acceleration v velocity along the y-axis G Grashof number D stretching factor j micro-inertia density V σ suction velocity N microrotation x distance along the sheet P r Prandtl number y distance normal to the sheet Sc Schmidt number s surface parameter p pressure b nonlinear (stretching) parameter q surface heat flux Greek η similarity variable γ micropolar parameter µ dynamic viscosity θ dimensionless temperature ν kinematic viscosity β coefficient of thermal expansion ν s microrotation/spin-gradient viscosity φ dimensionless mass transfer ρ density of the fluid λ suction parameter κ coupling coefficient ψ stream function

Subscripts
w surface conditions ∞ conditions far away from the surface

Introduction
In numerous industrial transport processes, convective heat and mass transfer takes place simultaneously.Phenomena involving stretching sheets feature widely in for example, aerospace component production metal casting [1].In such processes metals or alloys are heated until molten, poured into a mould or die, and liquid metal is subsequently stretched to achieve the desired product.When the super heated melt issues from the dies it loses heat and contract as it cools, a stage in metallurgical processing referred to as liquid state contraction.With further cooling and the loss of latent heat of fusion, the atoms of the metallic alloy lose energy and are bound tightly together in a regular structure.The mechanical properties of the final product depend to a great extent on the heat and mass transfer phenomena, the cooling rate, surface mass transfer rate etc.Much numerical research has been conducted in metal sheet flows including studies by Lait et al. [2], Goldschmit et al. who examined viscoplastic metal flows [3], Goldschmit [4] who provides a finite element methodology for general metal flow forming, and more recently by Cavaliere et al. [5].Some other important industrial applications of stretching sheet transport phenomena are the extrusion of a polymeric sheet from a die or the drawing of plastic films [6].During the manufacture of these sheets, the melt issues from a slit and is subsequently stretched to achieve the desired thickness.The mechanical properties of the final product are strongly influenced by the stretching rate and on the rate of cooling, both parameters which can be controlled by engineers.Polyester being a flexible material can tolerate filament surface stretching during the course of ejection; therefore the surface velocity deviates from being uniform.A linear stretching rate is all the more desirable, however in practice the stretching rate is observed to change during the process and is therefore generally nonlinear.Such extensional polymeric flows have attracted considerable attention in mathematical and experimental rheology since the 1960s.Ballman [7] considered the stretching flow of a Polystyrene melt.Ziabicki [8] provided a classical treatise on the physical mechanisms of polymer fiber flow dynamics in extensional sheet and other industrial flows.Denn and Marrucci [9] considered the stretching flow of viscoelastic suspensions with applications in plastics manufacturing.Denson [10] studied the stretching regime in polymeric processing.Newtonian stretching flows have also been studied in detail.Sakiadis [11] initiated the study of boundary layer flows over a continuous solid surface moving with constant speed for viscous fluids.Due to entrainment of ambient fluid, the situation represents a different class of boundary layer problem, which has a solution substantially different from that of boundary layer flow over a semi-infinite flat plate.For a Navier-Stokes (viscous) fluid, the heat and mass transfer on a stretching sheet with suction or blowing was investigated by Gupta and Gupta [12].Various scenarios in thermoconvective heat and mass transfer for stretching flows were subsequently discussed by many researchers.Surma Devi et al. [13] reported on numerical (finite difference solutions) for the transient three-dimensional boundary layer flow caused by a stretching surface.Chen and Char [14] studied the convection flow past a stretching surface with transpiration (wall mass flux) effects.Kumari et al. [15] modeled the hydromagnetic convection from a stretching surface with prescribed wall temperature/heat flux surface conditions using the Keller-Box numerical method.Surma Devi et al. [16] considered the momentum and heat transfer on a stretching sheet, presenting finite difference computations.Takhar and Nath [17] studied the three-dimensional flow due to a stretching flat surface with transient effects using shooting methods.They [18] extended this analysis to consider for the first time, rotational and magnetic body force effects on transient flow over a stretching surface.Takhar et al. [19] more recently analyzed the flow dynamics and species mass transfer in a stretching sheet with chemical reaction and magnetic retardation effects, using the Blottner difference scheme.
All the above investigations were restricted to Newtonian or non-Newtonian (viscoelastic or power-law type) fluids.In many environmental and industrial flows the classical theory of Newtonian fluids is unable to explain the microfluid mechanical characteristics observed.Micropolar fluids are fluids with microstructure belonging to a class of complex fluids with nonsymmetrical stress tensor referred to as micromorphic fluids.Physically they represent many industrially important liquids consisting of randomlyoriented particles suspended in a viscous medium.The classical theories of continuum mechanics are inadequate to explicate the microscopic manifestations of such complex hydrodynamic behaviour.Eringen [20] presented the earliest formulation of a general theory of fluid microcontinua taking into account the inertial characteristics of the substructure particles, which are allowed to sustain rotation and couple stresses.Later Eringen [21] generalized the theory to incorporate thermal effects in the so-called thermomicropolar fluid.The theory of micropolar fluids and its extension, the thermomicropolar fluid constitute suitable non-Newtonian hydrodynamic and thermo-hydrodynamic models which can simulate the flow dynamics of colloidal fluids, liquid crystals, polymeric suspensions, haemotological fluids etc.Many numerical studies of micropolar heat and mass transfer have been communicated in the literature.Hassanien and Gorla [22] investigated the heat transfer to a micropolar fluid from a nonisothermal stretching sheet with suction and blowing.Flow over a porous stretching sheet with strong suction or injection was examined by Kelson and Farell [23].Physically these investigations were confined to the linearly stretching sheet case.In certain polymeric processes nonlinear stretching effects are important.In metallurgical process the nonlinear velocity of a metal stream has to incorporated into models to provide a realistic computation.High velocity streams are erosive in nature and likely to dislodge particles of mould material.They are also likely to cause penetration of the metal in to mould material when they impinge upon it.On the other hand low velocity streams are more susceptible to aspiration and less likely to completely fill the mould cavity.Industrially therefore a moderate velocity is chosen to yield the desired characteristics.Drawing of plastic films and artificial fibers, are other applications of nonlinear stretching flows.
The purpose of the present investigation is therefore to study the coupled fluid flow, heat and mass transfer phenomena over a stretching sheet with nonlinear velocity, for micropolar fluids.Such a study has not been reported earlier in the literature and is important in non-Newtonian materials processing.It constitutes the first study to the author's knowledge of free convective flow with heat and mass transfer in micropolar nonlinear stretching flow regimes.The governing system of conservation equations is reduced to a coupled, multiple-degree system of nonlinear differential equations, which are solved by using the finite element method and also a finite difference method.Graphs are plotted for dimensionless translational velocity, micro-rotation (angular velocity), temperature and mass transfer function for various values of the nonlinear stretching parameter, Grashof number and Schmidt number.We also compute dimensionless wall shear stress function and heat transfer rate, demonstrating excellent agreement between the two numerical schemes employed.

Mathematical model
Consider an isothermal, steady, laminar, incompressible micropolar fluid flowing over a surface coinciding with the plate y = 0, the flow being confined in the region y > 0. Two equal and opposite forces are introduced along the x-axis so that the surface is stretched keeping the origin.The physical regime is illustrated in Fig. 1.The governing conservation equations, assuming thermal equilibrium can be cast as follows: where C A is the concentration of the microconstituent present in the flow.The thermal diffusivity k f in equation ( 4) and molecular diffusivity k g in equation ( 5) are assumed constant.Using the boundary layer concept, the effect of ∂N ∂x in the momentum equation is negligible in comparison to the effect of ∂N ∂y .It is similar to the fact that momentum equation for the velocity v in the y direction has been ignored.Similar considerations have been employed successfully by Kelson and Farell [23] and Bhargava et al. [24].The corresponding boundary conditions are given by Here C Aw > C A∞ , the transfer of species is occurring due to convection from the surface to the full stream fluid.A linear relationship between microrotation function N and surface shear ( ∂u ∂y ) is chosen for investigating the effect of different surface conditions for the microrotation.Here s is the boundary parameter and varies from 0 to 1.The first boundary condition (s = 0) is a generalization of the no-slip condition, which requires that the fluid particles closest to the solid boundary stick to it neither translating nor rotating.The second boundary condition i.e. microrotation is equal to the fluid vorticity at the boundary (s = 0), implies that in the neighborhood of a rigid boundary, the effect of microstructure is negligible since the suspended particles can not get closer to the boundary than their radius.Hence in the neighborhood of the fluid shear and therefore the gyration vector must be equal to fluid vorticity.
Introducing the dimensionless functions, f (η) and g(η), such that the continuity equation is automatically satisfied and choosing the similarity transformation as given by Vajravelu [25].
where prime denotes differentiation with respect to η.The governing equations ( 2) to ( 5) are then reduced to the following set of nonlinear, coupled ordinary differential equations.
The corresponding boundary conditions ( 6) reduce to: where is the local Grashof number i.e. convective parameter.b = 1 correspond to the linear stretching sheet.The parameters A and G correspond to local effects i.e pertaining to specific values of x.Similar studies were made by Kelson and Desseaux [26] and have been adopted in the present analysis.Equation ( 9) is therefore also an ordinary differential equation.P r = ρcpν k f is the Prandtl number and Sc = ρcpν kg is the Schmidt number.The shear stress at the surface of the sheet is defined as: The wall heat flux is computed using the following expression: The set of equations ( 8) to (11) are highly nonlinear and therefore the system cannot be solved analytically.The finite element method and also a finite difference technique have been used to solve the equations under the prescribed transformed boundary conditions (12a), (12b); both methods will now be described.

Numerical methods of solution 4.1 The finite element method
Finite element method is widely used for solving boundary value problems.The basic concept is that the whole domain is divided in to smaller elements of finite dimensions called "Finite Elements".Thereafter the domain is considered as an assemblage of these elements connected at a finite number of joint called "Nodes".The concept of discretization is adopted here.Other features of the method include seeking continuous polynomial, approximations of the solution over each element in term of nodal values, and assembly of element equations by imposing the inter-element continuity of the solution and balance of the interelement forces.The method entails the following steps: 1. Division of the domain in to linear elements, called the finite element mesh.
2. Generation of the element equations using variational formulations.
3. Assembly of the element equations as obtained in steps (2).
4. Imposition of the boundary conditions to the equations obtained in (3).

Solution of the assembled algebraic equations.
The assembled equations can be solved by any of the numerical technique viz.Gaussian elimination.An important consideration is that of shape functions which are employed to approximate actual functions.For one-dimensional and two-dimensional problems, the shape functions can be linear/quadratic and higher order.However the suitability of the shape functions varies from problem to problem.Due to the simple and efficient use in computations, linear as well quadratic shape functions are used in the present problem.However it is observed that the results do not vary considerably indicating that both elements provide approximately the same accuracy.The comparison for both types of shape functions is given in the Table 1.In our computations, the shape functions for a typical element (η e , η e+1 ) are taken as: linear element , η e ≤ η ≤ η e+1 ; quadratic element The general details of the steps employed in finite element analysis can be found in [27] and are summarized in Fig. 2 below.To solve the differential equations ( 8) to (11) with boundary condition (12), we assume: The system then reduces to: The corresponding boundary conditions (12a), (12b) reduce to: The whole domain is divided into eighty two-noded line elements, over each of the element, finite element equations are derived.
Variational formulation.The variational form associated with equations ( 15) to ( 19) over a typical linear element is given by ηe+1 ηe where w 1 , w 2 , w 3 , w 4 and w 5 are arbitrary test function and may be viewed as the variation in f, h, g, θ and φ respectively.All functions satisfy all homogeneous boundary conditions, as per theoretical requirements.

Finite element formulation.
As the domain is defined into two-noded elements, hence the appropriate finite element approximation is assumed as where w i = ξ 1 for the first node and w i = ξ 2 for the second node with i = 1, 2, 3, 4 and 5.Here ξ j are the shape functions for the line element (η e , η e+1 ) and are taken as: where η e ≤ η ≤ η e+1 .The finite element model of the equations ( 21) to (25) for the typical element is given by: Here each [K mn ] is of the order 2 × 2 and [r m ] (m, n = 1, 2, 3, 4, 5) is of 2 × 1.These matrices are defined as: where f = 2 j=1 f j ξ j , h = 2 j=1 h j ξ j .. The systems of equations after assembly of the elements are nonlinear therefore an iterative scheme is used to solve it.The system is linearized by incorporating the known function f and h.The whole domain is divided in to a set of 80 line elements.Each element matrix is of the order 10 × 10.Thus after assembly of all the elements equations we obtained a matrix of order 405 × 405.For the computational purpose η = ∞ has been fixed as 8.If η = ∞ is taken to be more than 8, all the unknown functions do not change up to the desired accuracy.

Finite difference method
For comparison purposes the same system of equations ( 15)-( 19), subject to boundary conditions (20) are solved numerically using the finite difference method.This method is used for solving ordinary as well as partial differential equations governing boundary value problem as well as initial value problem.This method can be explained briefly the following Fig. 3: By using the central difference formulae, the set of equations ( 15)- (19), can be written as: where h e is the step length.Since the above equations are non-linear and coupled hence they cannot be solved exactly.Therefore an iterative scheme is required to be used.Writing down the equations in the form: where each l i is the function of the variable f i , h i , g i , θ i and x i is any of the variable f i , h i , g i , θ i .Similarly equations are formulated for each variable of the equations ( 30)-(34).Commencing with the initial guess values, new iterate values are obtained.This process continues until the absolute error |x i − x i−1 | is less than the accuracy required.The condition of convergence of the scheme has been already checked before implementing the iterative scheme.Following equation ( 35), the equations ( 30)-(34) can be written as follows: The boundary conditions are presented as: The system of equations ( 36) to (40) with the boundary conditions (41) has been solved iteratively and the results obtained are compared with those obtained by FEM.

Results and discussion
The variation of the skin friction and heat transfer with respect to convective parameter G, surface parameter s, Prandtl number P r, and nonlinear parameter b are depicted in Table 2.It can be seen that the skin friction coefficient −f (0) ′′ , increases with increase in s, nonlinear parameter and Prandtl number while it decreases with buoyancy parameter (G).The rate of heat transfer, −θ(0) ′ , increases with higher values of b and Prandtl number and it decreases with an increase in convective parameter G.   Comparison between the finite element and finite difference solutions is illustrated in Table 3, where for s = 0.5, P r = 7, b = 0.5, Sc = 1.0,G = 0.5 we have compared profiles of h, g and θ with η coordinate.Excellent correlation is demonstrated between the two numerical methods.We observe that h (dimensionless translational velocity), g (dimensionless micro-rotation) and θ (dimensionless temperature) all decrease continuously from a peak value of unity at η = 0 to a minimum value at η = 8.In addition we have computed these profiles using both linear and quadratic elements with the finite element program, again for arbitrary values of the thermophysical parameters and observe very little difference in the computations.
In the graphs provided, the dimensionless velocity, microrotation, temperature and mass transfer functions are computed for fixed value of Prandtl number P r, material parameter α 1 , surface parameter s, and physical micropolar parameters, A and A 1 , which are taken as 7.0, 1.0, 0.5, 1.0 and 1.0 respectively, while the effect of other important parameters namely, nonlinear parameter, b, Grashof number, G and Schmidt number, Sc, have been studied explicitly.These functions are shown in Figs. 4, 5.
Table 3.Comparison of the results (FEM -finite element Vs, FDM -finite difference).s = 0.5, P r = 7., b = 0.5, Sc = 1.0,G = 0.This indicates that with strong nonlinear stretching micropolar fluids achieve a decrease in temperature compared to Newtonian fluids, which may be beneficial in temperature control of polymer stretching processes.Generally as b rises the temperature increases with an increase in the nonlinear stretching parameter b.Fig. 4(d) shows the effect of the stretching parameter, b, on mass transfer function, φ, for the case of a micropolar fluid (R = 0.5).Mass transfer function, φ, clearly increases markedly with a rise in b; as expected fast stretching therefore enhances mass transfer.All profiles descend smoothly from unity at the wall to zero in the free stream (far from the wall).Furthermore the variation in all above functions becomes less pronounced for large b values.This is due to the fact that the coefficient 2b b+1 , 1 b+1 in the differential equations (8) and 3b−1 b+1 in (9) approaches to 2, 0 and 3 respectively as b approaches infinity, which is true for real flows as the stretching parameter cannot be too large i.e. must have a limiting value.Fig. 5 illustrates the variation of dimensionless translational velocity, microrotation, temperature and mass transfer functions with η for various Grashof numbers, G. Fig. 5(a) demonstrates the variation of the velocity distribution with the free convective parameter, G.It is observed that velocity continuously increases with an increase in the Grashof number G, which implies an increase in the boundary layer thickness.Velocity is observed to be a maximum near the boundary and decreases far away from the boundary.For moderate value of G the velocity profile changes its nature; it exhibits a steep behaviour and decreases as Grashof number decreases.Increasing buoyancy (i.e.G values) therefore enhance the velocity values i.e. accelerate the flow.Fig. 5(b) shows the variation of the microrotation distribution, g, which continuously decreases with an increase in the Grashof number, G.For large values of the convective parameter, G, near the boundary, microrotation is negative (i.e. for G = 20), whereas away from the boundary it becomes positive and finally descends to zero.Thus large convection effects produce a reverse rotation only near the boundary.All microrotation profiles for G = 0.5, 1, 10 and 20, converge at η ∼ 3 and tend to zero.The temperature distribution, θ, with η, is shown in Fig. 5(c).θ increases continuously with an increase in the value of the Grashof number, G.The convective parameter has thus an important role in controlling the temperature.As in many metallurgical processes [1], temperature rises to a great value (during intermediate stages), and there is a need of maintain appropriate temperatures, which can be achieved using convective (buoyancy) effects.Lower Grashof numbers therefore depress temperatures throughout the flow regime.In Fig. 5(d) we observe that the mass transfer function, φ, is increased slightly with a rise in convection parameter, G. Hence the largest values of φ correspond to G = 20, and these decrease as G falls to 10, 1 and 0.5.Mass transfer therefore may also be inhibited by reducing the Grashof number, G, in practical applications.Finally in Fig. 6 we have illustrated the profile of mass transfer function, φ, with η coordinate for various Schmidt numbers, Sc.All profiles descend smoothly from unity at the wall (η = 0) to zero in the free stream (η = 8).A rise in Sc from 0.5 through 1, 2, 3 to 5 induces a considerable reduction in φ values indicating that mass transfer is reduced in the micropolar fluid with higher Sc values.For higher Sc values, the profiles also descend much faster to zero; this descent becomes more gradual as Sc decreases in value.

Conclusions
1.The dimensionless translational velocity, f ′ , for both micropolar and Newtonian fluids decreases with an increase in nonlinear parameter b, although values for micropolar fluids are consistently higher.For micropolar fluids (using in our computations, R = 0.5), f ′ increases with an increase in the convective parameter G, while it.
2. Dimensionless micro-rotation (angular velocity), g, generally increases with a rise in nonlinear parameter b (in particular in the vicinity of the wall) and decreases with a rise in convective parameter, G (again most notably in the vicinity of the wall).With considerable distance from the wall (η > 3), G however has negligible effect on micro rotation profiles.
3. Convective parameter, G i.e.Grashof number, can be used effectively for controlling the temperature field.5. Dimensionless mass transfer function, φ, for micropolar fluids, increases generally with an increase in nonlinear parameter, b and also with an increase in convective parameter, G; however φ decreases with a rise in Schmidt number.
6. Skin friction numerically increases with an increase in surface parameter s, Prandtl number P r and nonlinear parameter b, while decreases with increase in convective parameter G, indicating that these parameters accelerate the flow regime.Thus drag forces can be reduced by an increase in the free convection parameter, G.
7. Heat transfer rate increases with an increase in convective parameter G, nonlinear stretching parameter, b, and also Prandtl number.The increase in heat transfer rate indicates a fast cooling of the plate.However a rise in surface parameter, s, is shown to decrease the heat transfer rate markedly.

Fig. 4
Fig.4illustrates the variation of dimensionless translational velocity, microrotation, temperature and mass transfer (concentration) functions with nonlinear parameter b.Profiles are depicted both for R = 0 and R = 0. R = 0 corresponds to the Navier-Stokes viscous fluid, which was considered by Vajravelu[25].Fig.4(a) demonstrates the variation of velocity with the parameter b, where velocity decreases as b increases.It is found that the translational velocity for a Newtonian fluid is less than that for micropolar fluids.The results for b = 1 correspond to a linear stretching sheet, for which the lowest translational velocities are observed.Fig. 4(b) shows the variation of the microrotation function, g, with b.It is clear from the figure that microrotation function first increases near the boundary as b increases, but away from the boundary a reverse pattern is observed.Fig. 4(c) depicts the variation of temperature function, θ, with b.Temperature for Newtonian flow is less than the micropolar fluids for b = 1 i.e. the linear stretching case.However for b = 5 the temperature values are almost identical for both Newtonian (R = 0) and micropolar (R = 0.5) cases.For the highest value of b (strong nonlinear stretching) we observe a reversal from the linear case (b = 1), wherein micropolar fluid temperature (R = 0.5) is less than Newtonian fluid temperature.This indicates that with strong nonlinear stretching micropolar fluids achieve a decrease in temperature compared to Newtonian fluids, which may be beneficial in temperature control of polymer stretching processes.Generally as b rises the temperature increases with an increase in the nonlinear stretching parameter b.Fig.4(d)shows the effect of the stretching parameter, b, on mass transfer function, φ, for the case of a micropolar fluid (R = 0.5).Mass transfer function, φ, clearly increases markedly with a rise in b; as expected fast stretching therefore enhances mass transfer.All profiles descend smoothly from unity at the wall to zero in the free stream (far from the wall).Furthermore the variation in all above functions becomes less pronounced for large b values.This is due to the fact that the coefficient 2b b+1 , 1 b+1 in the differential equations (8) and 3b−1 b+1 in (9) approaches to 2, 0 and 3 respectively as b approaches infinity, which is true for real flows

4 .
Dimensionless temperature, θ for both micropolar and Newtonian fluids, increases with an increase in nonlinear parameter, b.For micropolar fluids our results indicate that θ increases with a rise in convective parameter, G.

Table 1 .
Comparison of results with linear as well as quadratic elements

Table 2 .
Table for skin friction −f ′′ (0) and the rate of heat transfer −θ ′ (0) with different value of surface parameter s, Grashof number G, Prandtl number P r and the nonlinear parameter b s