Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 52, 2, pp. 557-569, Warsaw 2014 NUMERICAL SIMULATION OF THE MICROPOLAR FLUID FLOW AND HEAT TRANSFER IN A CHANNEL WITH A SHRINKING AND A STATIONARY WALL Kashif Ali, Muhammad Ashraf Centre for Advanced Studies in Pure and Applied Mathematics, Bahauddin Zakariya University, Multan, Pakistan e-mail: kashifali 381@yahoo.com We present the numerical study of hydromagnetic (MHD) flow and heat transfer characte- ristics of a viscous incompressible electrically conductingmicropolar fluid in a channel with one wall shrinking and the other at rest in the presence of a transverse applied magnetic field. Different from the classical shooting methodology, we employ a combination of a di- rect andan iterativemethod (SORwith optimal relaxationparameter) for solving the sparse systems of linear algebraic equations arising from the FDdiscretization of the linearized self similar nonlinear ODEs. Effects of some physical parameters on the flow and heat transfer are discussed and presented through tables and graphs. The present investigation may be beneficial to the flow and thermal control of polymeric processing. Keywords: channel flow, shrinking/stationarywalls, couple stress, quasi-linearization, Joule heating, viscous dissipation 1. Introduction Flows through channels have attracted the research community due to their applications in the fields of binary gas diffusion, microfluidic devices, surface sublimation, ablation cooling, filtration, grain regression (as in the case of combustion in rocket motors) and the modeling of air circulation in the respiratory system (see Ashraf et al., 2009b). Laminar air-flow systems have been used by the aerospace industry to control particulate contamination. Furthermore, the laminar flow cabinets have been used in the maintenance of negative pressure and in the adjustmentof the fans to exhaustmore air.Therefore, theNavier Stokes equations,whichare the governing equations for these problems, have attracted the interest of the researchers. An exact solution of Navier Stokes equations for motion of an incompressible viscous fluid in a channel with different pressure gradients was presented by Suton andBarto (2008). Their solutions were helpful in verifying and validating computational models of complex unsteady flows to guide the design of fuel injectors and controlled experiments. Numerical simulation of flow through microchannels with design roughness was presented byRawool et al. (2006). Recently, Hajipour and Dehkordi (2012) studied the transient behavior of fluid flow and heat transfer in vertical channels partially filled with a porousmedium including the effects of inertial term and viscous dissipation. All the studies cited above are limited to classical Newtonian fluids. There are many fluids which are important fromthe industrial point of view, anddisplaynon-Newtonian behavior.Due to the complexity of such fluids, several models have been proposed but the micropolar model is the most prominent one. Hoyt and Fabula (1964) experimentally predicted that the fluids which could not be characterized by Newtonian relationships indicated significant reduction in shear stress near a rigidbodyand couldbewell explainedby themicropolarmodel introducedby Eringen (1964).MHDflowandheat transfer of anon-Newtonianfluidover a stretching/shrinking surface is important due to its engineering and industrial applications. For example, in the 558 K. Ali, M. Ashraf extrusion of a polymer sheet fromadye, the sheet is sometimes stretched,whereas the properties of the end product depend considerably on the rate of cooling. The micropolar fluid is a hot area of research and thereforemany investigators have studied the related flowand heat transfer problems. For example, Shangjun et al. (2006), Kelson et al. (2003), Naccache andSouza (2011), Ashraf et al. (2009a,b, 2011) and Ashraf and Batool (2013). In the present paper, the hydromagnetic flow and heat transfer characteristics of a visco- us electrically conducting incompressible micropolar fluid in a channel with a shrinking and a stationary wall are presented. Viscous dissipation and Joule heating effects have also been taken into account. A similarity transformation is used to convert the governing Navier-Stokes equations into a set of nonlinear ODEs, which are numerically solved using a computational procedure based on the quasi-linearization and finite difference discretization. The effects of the governing parameters on the flow and heat transfer aspects of the problem are discussed. 2. Mathematical formulation Weconsider two dimensional hydromagnetic steady laminar viscous flow and heat transfer of an incompressiblemicropolar fluid in a parallel plate channel with the lower wall shrinking and the upper at rest, in the presence of a transverse appliedmagnetic field. The inducedmagnetic field is assumed to be negligible as comparedwith the imposed field. Themagnetic Reynolds number which is the ratio of the product of characteristic length and fluid velocity to the magnetic diffusivity is assumed to be small. It is used to compare the transport of magnetic lines of the force in a conducting fluid, with the leakage of such lines from the fluid. For a small magnetic Reynolds number, the magnetic field will tend to relax towards a purely diffusive state. It is also assumed that there is no applied polarization and hence no electric field. The body couple is assumed to be absent. The two walls of the channel are located at y=−a and y= a, where 2a is the channel width. We, therefore, may express the velocity v and the microrotationυ, respectively, in the following form v=(u(x,y),v(x,y),0) υ=(0,0,φ(x,y)) (2.1) where φ is the component of the microrotation normal to the xy-plane. Taking the viscous dissipation and Joule heating effects into consideration, the governing equations of the flow and heat transfer for the problem are ∂u ∂x + ∂v ∂y =0 u ∂u ∂x +v ∂u ∂y = −1 ρ ∂p ∂x + µ+κ ρ ∇2u+ κ ρ ∂φ ∂y − σB20 ρ u u ∂v ∂x +v ∂v ∂y = −1 ρ ∂p ∂y + µ+κ ρ ∇2v− κ ρ ∂φ ∂x ρj ( u ∂φ ∂x +v ∂φ ∂y ) = γ∇2φ+κ (∂v ∂x − ∂u ∂y ) −2κφ ρcp ( u ∂T ∂x +v ∂T ∂y ) = k0 ∂2T ∂y2 +µ (∂u ∂y )2 +σB20u 2 (2.2) where ρ is the density, p is the pressure, µ is the dynamic viscosity, κ is the vortex viscosity, σ is the electrical conductivity, B0 is the strength of the magnetic field, j is the micro-inertia density, γ is the spin gradient viscosity, cp is the specific heat at constant pressure, κ0 is the thermal conductivity and T is the temperature of the fluid. The boundary conditions for the velocity, microrotation and temperature fields for the present problem are Numerical simulation of the micropolar fluid flow and heat transfer... 559 u(x,−a)=−bx u(x,a)= 0 v(x,±a) = 0 ϕ(x,±a) = 0 T(x,−a)=T1 T(x,a) =T2 (2.3) where b > 0 is the shrinking rate of the channel wall, and T1 and T2 (with T1 > T2) are fixed temperatures of the lower and upper channel walls, respectively. We define the following similarity variables to convert governing partial differential Eqs. (2.2) into ordinary differential equations η= y a u= bxf ′(η) v=−abf(η) ϕ=− b a xg(η) θ(η)= T −T2 T1−T2 (2.4) The above proposed velocity field is compatible with continuity Eq. (2.1) and, therefore, repre- sents the possible fluidmotion. After eliminating the pressure term fromEqs. (2.2)2 and (2.2)3, we introduce the above similarity transformation to the resulting equation to get (1+C1)f iν −C1g ′′ =Re(f ′f ′′−ff ′′′)+Mf ′′ (2.5) whereas Eqs. (2.2)4 and (2.2)5, in view of Eq. (2.4), are C3g ′′+C1C2(f ′′−2g) = f ′g−fg′ θ′′+PrRefθ′+EcPr(f ′′ 2 +Mf ′ 2 )= 0 (2.6) Here C1 = κ/µ is the vortex viscosity parameter, C2 =µ/(ρjb) is the microinertia density pa- rameter, C3 = γ/(ρja 2b) is the spin gradient viscosity parameter, Re= ρa2b/µ is the Reynolds number, Ec = (bx)2/[cp(T1 −T2)] is the Eckert number, Pr = µcp/κ0 is the Prandtl number and M = a2σB20/µ is themagnetic parameter. Boundary conditions (2.3), in view of Eq. (2.4), in dimensionless form are reduced to f(±1)= 0 f ′(−1)=−1 f ′(1)= 0 g(±1)= 0 θ(−1)= 1 θ(1)= 0 (2.7) 3. Numerical solution A usual approach is to write the nonlinear ODEs in form of a first order initial value problem as follows: Setting: f ′ = p, f ′′ = q, f ′′′ = r, g′ = s, θ′ = t in Eqs. (2.5) and (2.6), we have f ′ = p, p′ = q q′ = r r′ = 1 1+C1 (Repq−Refr+Mq)+ C1 (1+C1)C3 (pg−fs−C1C2q+2C1C2g) g′ = s s′ = 1 C3 (pg−fs−C1C2q+2C1C2g) θ′ = t t′ =−PrRetf −EcPr(q2+Mp2) (3.1) with the following required boundary conditions f(−1)= 0 p(−1)=−1 g(−1)= 0 θ(−1)= 1 q(−1)=α1 r(−1)=α2 s(−1)=α3 t(−1)=α4 (3.2) Here, α1,α2,α3,α4 (that is f ′′(−1), f ′′′(−1), g′(−1), θ′(−1)) are theunknown initial conditions. Therefore, a shooting methodology is incorporated to solve the above system, which may be a 560 K. Ali, M. Ashraf combination of the Runge Kuttamethod (to solve the 1st order ODEs) and a four dimensional zero finding algorithm (to find the missing conditions). It is worthy to note that the missing initial conditions are computed so that the solution satisfies the boundary conditions f(1)= 0, p(1)= 0, g(1)= 0, θ(1)= 0, of the original boundary value problem. A serious shortcoming of the shooting is the blowing up of the solution, before the initial value problem is completely integrated, and this happens quite often even with very accurate guesses for the initial conditions. This phenomenon is due to the instability of the differential equations and also because of the inherent strong dependence of the solution on the initial conditions of the problem.On the other hand, a finite differencemethod (FDM) does not suffer from this short-coming, and does have a chance as it tends to keep a firm hold on the entire solution at once. Due to these features, it is desirable to incorporate theFDMin the computational procedure. Therefore, in our earlier works (Ashraf et al., 2009a,b, 2011; Ashraf and Batool, 2013), we followed an approach different from the usual shooting method and involving the FDM, which we describe as follows: An inspection of Eq. (2.5) reveals that it can be readily integrated, and becomes (1+C1)f ′′′−C1g ′ −Re(f ′ 2 −ff ′′)=β (3.3) where β is the constant of integration to be determined. We consider this third order equation as the following system of coupled 1st and 2nd order ODEs f ′ = p (1+C1)p ′′−C1g ′−Re(p2−fp′)=β (3.4) We solve these equations along with Eqs. (2.6), subject to the boundary conditions f(±1)= 0 p(−1)=−1 p(1)= 0 g(±1)= 0 θ(−1)= 1 θ(1)= 0 (3.5) For the numerical solution of the present problem, we first discretize Eqs. (2.6) and (3.4)2 at a typical grid point η= ηi by employing central difference approximations for the derivatives. The resultant algebraic system is solved iteratively by the SORmethod, subject to the appropriate boundary conditions given in Eq. (3.5). On the other hand, Eq. (3.4)1 is integrated numerically by employing the 4th order Runge Kuttamethod, after every SOR iteration. It is important tomention that theboundarycondition fn =0 is used in finding the constant of integration β which is the only unknown as compared to the four missing initial conditions in the shooting approach. In our earlier work, we used the trial and errormethod to find β but some one dimensional zero finding algorithms may also be employed for the purpose. However, the sensitivity of β with respect to the governing parameters makes it difficult to find it, and some manual effort is always required in every simulation to obtain its desired value. That is why we need some alternative approach which does not require finding any unknown and is entirely based on the FDM. In this paper, we discuss an alternative approach based on the quasi-linearization of nonlinear ODEs. 3.1. Quasi-linearization method In the quasi-linearization, we construct sequences of vectors {f(k)}, {g(k)}, {θ(k)} which converge to the numerical solutions of Eqs. (2.5) and (2.6), respectively. To construct {f(k)}, we linearize Eq. (2.5), by retaining only the first order terms, as follows: Numerical simulation of the micropolar fluid flow and heat transfer... 561 We set G(f,f ′,f ′′,f ′′′,fiv)≡ (1+C1)f iv+Reff ′′′−Ref ′f ′′−Mf ′′−C1g ′′ G(f(k),f(k) ′ f(k) ′′ ,f(k) ′′′ ,f(k) iv )+(f(k+1)−f(k)) ∂G ∂f(k) +(f(k+1) ′ −f(k) ′ ) ∂G ∂f(k) ′ +(f(k+1) ′′ −f(k) ′′ ) ∂G ∂f(k) ′′ +(f(k+1) ′′′ −f(k) ′′′ ) ∂G ∂f(k) ′′′ +(f(k+1) iv −f(k) iv ) ∂G ∂f(k) iv =0 (1+C1)f (k+1)iv+Ref(k)f(k+1) ′′′ − (Ref(k) ′ +M)f(k+1) ′′ −Ref(k) ′′ f(k+1) ′ +Ref(k) ′′′ f(k+1) =Re(−f(k) ′ f(k) ′′ +f(k)f(k) ′′′ )+C1g (k)′′ (3.6) Now Eqs. (3.6) gives a system of linear differential equations with fk being the numerical solution vector of the k-th equation. To solve the ODEs, we replace the derivatives with their central difference approximations, giving rise to the sequence {f(k)}, generated by the following linear system Af(k+1) =B with A=A(f(k)) and B=B(f(k),g(k)) (3.7) where n is the number of grid points. Thematrices An×n and Bn×1 are initialized as follows A1,1 =1 B1 =0 A2,1 =−4(1+C1)+hRef1−Mh 2 A2,2 =7(1+C1)−hRef2+h 2Reα+ hRe 2 f4+2Mh 2 A2,3 =−4(1+C1)−hRef3−Mh 2 A2,4 =1+C1+ hRe 2 f2 B2 =h 2C1(g3 −2g2+g1)+ hRe 2 f2(−f2+2f1−2f3+f4) − hRe 2 (f3−f1)(f3−2f2+f1)+2hα(1+C1) (3.8) and for 2 1 and µ (t) ω < 1. The iterative process is continued with ω = ωopt until convergence is achieved. TOLωopt is chosen to lie in the range 10 −5 ¬ TOLωopt ¬ 10 −3 depending upon the density of the grid points. For a coarse grid, larger values result in rapid convergence while for a fine grid, smaller values may be chosen for better estimation of ωopt. We may improve the order of accuracy of the solution by solving the problem again on the grid with step h/2 and h/4, and then using the Richardson extrapolation which have been carried out at not only the common grid points but also at the skipped points, by following Roache and Knupp (1993). Any extrapolation scheme (Deuflhard, 1983) may be used for this purpose but we have used the polynomial extrapolation. Moreover, a good initial guess for the solution on the finer grid can be obtained by injecting the previous coarse grid solution. Many operators may be found for this purpose but we have used the following linear operator which is simpler and easier to implement W2i−1 ←U c i W2i ← 1 2 (Uci +U c i+1) (3.16) where Uc represents the coarse grid solution and W is the required initial guess. 4. Results and discussion The physical quantities of our interest are the shear stresses, the couple stresses and the heat transfer rates at the channel walls which are, respectively, proportional to the values of f ′′, g′, θ′ at the two walls. The parameters of the problem are the Reynolds number Re, the magnetic parameter M, the Prandtl number Pr, the Eckert number Ec and the micropolar parameters C1, C2 and C3. These parameters are all dimensionless groups of material and flow properties, and/or geometric dimensions of the domain. The traditional way (e.g., Lok et al., 2007) of stu- dying flow and thermal characteristics of the fluid dynamics problems is to specify the values of these dimensionless groups rather than specifying the particular fluid properties and the doma- in dimensions. Obviously, the results obtained in this way are applicable to the flow problems with particular values of material properties and the dimensions of the domain, falling in the ranges considered in the studies. In the present work we have computed our results for the fol- lowing values of the governing parameters: M =0,10,20,30,40,50, Re= 1,5,10,20,30,40,50, Ec=0,1,2,3,4 and Pr=0.1,0.5,1,2,2.5. We shall study the effects of the parameters on f ′′(±1), g′(±1) and θ′(±1) as well as on the velocity profiles f(η), f ′(η), the microrotation profile g(η) and the heat profile θ(η). All the cases of themicropolar parameters in the presentwork are shown inTable 1, whereas thefirst case corresponds to theNewtonianfluid.Table 2 shows the convergence of our numerical results as the step-size decreases, which gives us confidence on our computational procedure. The magnetic field exerts a force, called the Lorentz force, which tends to drag the fluid towards the shrinking wall, which not only results in increasing the shear stress at the wall but also causes greater spinning of the micro fluid particles, and hence increases the couple stress as well (clear from Table 3). Moreover, the frictional force raises the temperature of the 564 K. Ali, M. Ashraf Table 1. Five cases of values of micropolar parameters C1,C2 and C3 Cases C1 C2 C3 1(Newtonian) 0 0 0 2 2 0.7 0.3 3 4 1.5 0.7 4 6 2.0 1.5 5 8 2.5 2.0 Table 2. Dimensionless temperatures θ(η) on three grid sizes and extrapolated values for C1 =4,C2 =1.5, C3 =0.7, Re=50,M =1, Ec=0.5 and Pr=0.5 θ(η) η 1st grid η 2nd grid 3rd grid Extrapolated (h=0.02) (h=0.01) (h=0.005) values −0.8 1.015900 1.015916 1.015920 1.015921 −0.6 1.028347 1.028372 1.028379 1.028381 −0.4 1.039211 1.039244 1.039252 1.039255 −0.2 1.048684 1.048721 1.048731 1.048734 0 1.056010 1.056130 1.056137 1.056140 0.2 1.056186 1.056160 1.056153 1.056151 0.4 1.018667 1.018485 1.018439 1.018424 0.6 0.856763 0.856441 0.856361 0.856334 0.8 0.497270 0.497052 0.496996 0.496978 Table 3.The effect of themagnetic field on shear and couple stresses and the heat transfer rate with C1 =3,C2 =0.8, C3 =1, Re=20, Ec=0.5 and Pr=0.1 M f ′′(−1) f ′′(1) g′(−1) g′(1) θ′(−1) θ′(1) 0 0.297547 −2.871344 1.750801 3.082036 0.080047 −6.683872 10 0.663511 −2.214353 2.239587 2.458181 2.738352 −7.638689 20 1.082883 −1.736506 2.711143 1.970446 5.237685 −8.024192 30 1.509880 −1.412842 3.125411 1.616686 7.576004 −8.148545 40 1.915830 −1.197936 3.472198 1.366667 9.766241 −8.181408 fluid which results in increasing the temperature difference between the channel walls and the fluid. Therefore, the heat transfer rate at the two walls, which is directly proportional to the temperature difference, also increases. Table 4 predicts that the shear and couple stresses as well as the heat transfer rate at the shrinking wall may decrease as the Reynolds number increases. The increased shrinking rate of the wall forces the fluid tomove rapidly away from the shrinkingwall, thus decreasing both the shear and couple stresses. Also, the fluid carries away heat (generated due to the Joule heating effect) from the region near the lowerwall, thus decreasing the temperature difference andhence the heat transfer rate at the shrinkingwall. On the other hand, the incoming of the heated fluid towards the upperwall increases not only the shear and couple stresses but also the heat transfer rate. It is clear from Table 5 that the micropolar structure of the fluid tends to reduce the shear stress at the shrinking wall without significantly altering the shear stress at the stationary wall. This is in accordance with the experimental prediction of Hoyt and Fabula (1964) that the micro fluid particles cause significant reduction in the shear stress near a moving rigid surface. Moreover, the particles also cause microrotation in the fluid which is responsible for the couple Numerical simulation of the micropolar fluid flow and heat transfer... 565 Table 4.The effect of the Reynolds number on shear and couple stresses and the heat transfer rate with C1 =3,C2 =0.8,C3 =1,M =5, Ec=0.5 and Pr=0.1 Re f ′′(−1) f ′′(1) g′(−1) g′(1) θ′(−1) θ′(1) 1 2.048361 −0.780919 3.566882 1.130933 7.995645 −3.260171 10 1.755744 −0.941748 3.330322 1.306377 5.654433 −4.520490 20 1.431412 −1.174269 3.045316 1.545727 4.325219 −5.211422 30 1.127477 −1.471094 2.752210 1.830289 3.562179 −5.952669 40 0.867088 −1.825513 2.476137 2.143552 3.072278 −6.776226 Table 5.The effect ofmicropolar coefficients on shear and couple stresses and the heat transfer rate with Re=5,M =20, Pr=0.5 and Ec=0.5 Case f ′′(−1) f ′′(1) g′(−1) g′(1) θ′(−1) θ′(1) 1 4.516739 −0.711734 0 0 0.759419 −0.738241 2 2.895274 −0.764749 2.720250 0.561054 0.839335 −0.900799 3 2.365865 −0.773087 3.499924 1.040502 0.910777 −0.979324 4 2.149622 −0.774630 3.576582 1.049530 0.953960 −1.029429 5 1.981817 −0.781593 3.871003 1.273793 0.986792 −1.061255 Table 6. The effect of the Eckert number on the heat transfer rate with C1 = 3, C2 = 0.8, C3 =1, Re=10, M =1 and Pr=0.1 Ec θ′(−1) θ′(1) 0 −0.393561 −0.581297 1 −0.235738 −0.662955 2 −0.077915 −0.744612 3 0.079908 −0.826269 4 0.237731 −0.907926 Table 7.The effect of thePrandtl number on heat transfer rate with C1 =3,C2 =0.8,C3 =1, Re=10,M =1 and Ec=3 Pr θ′(−1) θ′(1) 0.1 0.079908 −0.826269 0.5 1.909219 −2.200202 1 3.321423 −3.787488 2 4.990100 −6.302292 2.5 5.615060 −7.350247 stress at the walls. Themicropolar structure also tends to increase the heat transfer rate at the two walls by raising the temperature distribution across the channel. On the other hand, it is evident fromTable 6 that the viscous dissipationmay cause thermal reversal near the shrinking wall only, which is supported by the Prandtl number as predicted in the Table 7. The streamlines for the present problem are shown in Fig. 1. On the other hand, Figs. 2 and 3 show that the external magnetic field reduces the two velocities (normal and streamwise) while raising the fluid temperature, whereas it raises themicrorotation profiles only in a smaller region near the shrinking wall. Moreover, it supports the thermal reversal near the shrinking wall, due to the viscous dissipation and Joule heating effects. Figures 4 and 5 show that the influence of the Reynolds number on the problem is almost opposite to that of the magnetic field, with temperature distribution near the upper wall being 566 K. Ali, M. Ashraf Fig. 1. Streamlines for the problem for C1 =4,C2 =1.5,C3 =0.7, R=1 and M =50 Fig. 2. Streamwise (a) and normal (b) velocity profiles for C1 =4,C2 =1.5,C3 =0.7, Re=50, Ec=1, Pr=2.5 and various M Fig. 3. Microrotation (a) and temperature (b) profiles for C1 =4,C2 =1.5,C3 =0.7, Re=50, Ec=1,Pr=2.5 and various M Fig. 4. Normal (a) and streamwise (b) velocity profiles for C1 =4,C2 =1.5,C3 =0.7,M =10, Ec=1, Pr= 2.5 and various Re Numerical simulation of the micropolar fluid flow and heat transfer... 567 Fig. 5. Microrotation (a) and temperatue (b) profiles for C1 =4,C2 =1.5,C3 =0.7,M =10, Ec=1, Pr= 2.5 and various Re Fig. 6. Normal (a) and streamwise (b) velocity profiles for M =20, Ec=1, Re=5, Pr= 2.5 and various cases of the micropolar parameters Fig. 7. Microrotation (a) and temperatue (b) profiles for M =20, Ec=1, Re=5, Pr=2.5 and various cases of the micropolar parameters Fig. 8. Temperature profiles for: (a) C1 =4,C2 =1.5,C3 =0.7, Re=50, Pr= 2.5,M =20 and various Ec, (b) C1 =4,C2 =1.5,C3 =0.7, Re=50, Ec=1,M =20 and various Pr 568 K. Ali, M. Ashraf the only exception. It is clear from Figs. 6 and 7 that the micropolar structure of the fluid increases the velocity, microrotation and temperature distributions across the channel. Finally, Fig. 8a shows that the viscous dissipation may cause thermal reversal near the shrinking wall by significantly raising the temperature profiles, which is supported by the Prandtl number as predicted by the Fig. 8b. 5. Conclusions We have numerically studied the problem of a viscous incompressible electrically conducting micropolar fluid in a channel with onewall shrinking and the other at rest. It has been observed that the externalmagnetic field tends to raise themicrorotation profiles only in a smaller region near the shrinking wall. Moreover, it tends to balance the influence of the shrinking channel wall on the flow and heat transfer characteristics of the problem. The combined effect of viscous dissipation and Joule heating may raise the temperature distribution to such an extent that thermal reversal may occur near the shrinking wall, which is supported by both the Prandtl number and themicropolar structure of the fluid. Acknowledgements The authors are extremely grateful to theHigher EducationCommission of Pakistan for the financial support to carry out this research. The authors are also very grateful to the learned reviewer for the useful comments improving the quality of the paper. References 1. AshrafM.,BatoolK., 2013,MHDflowandheat transfer of amicropolarfluid over a stretchable disk, Journal of Theoretical and Applied Mechanics, 51, 25-38 2. AshrafM.,KamalM.A., SyedK.S., 2009a,Numerical simulationof amicropolarfluidbetween a porousdisk andanon-porousdisk,Journal of AppliedMathematics andModelling, 33, 1933-1943 3. Ashraf M., Kamal M.A., Syed K.S., 2009b, Numerical study of asymmetric laminar flow of micropolar fluids in a porous channel,Computers and Fluids, 38, 1895-1902 4. AshrafM.,KamalM.A., SyedK.S., 2011,Numerical simulation of flowofmicropolar fluids in a channel with a porous wall, International Journal for Numerical Methods in Fluids, 66, 906-918 5. Deuflhard P., 1983,Order and step-size control in extrapolationmethods,Numerical Mathema- tics, 41, 399-422 6. Eringen A.C., 1964, Simple micropolar fluids, International Journal of Engineering Science, 2, 205-217 7. HajipourM.,DehkordiA.M., 2012,Transientbehaviorof fluidflowandheat transfer in vertical channel partially filledwith porousmedium:Effects of inertial termandviscous dissipation,Energy Conversion and Management, 61, 1-7 8. Hoyt J.W., Fabula A.G., 1964, The effect of additives on fluid friction, US Naval Ordinance Test Station Report 9. Kelson N.A., Desseaux A., Farrell T.W., 2003, Micropolar flow in a porous channel with high mass transfer,ANZIAM, 44, 479-495 10. LokY.Y., IoanP., ChamkhaA.J., 2007,Non-orthogonal stagnation point flow of amicropolar fluid, International Journal of Engineering Science, 45, 173-184 11. Naccache M.F., SouzaP.R., 2011,Heat transfer to non-Newtonian fluids in laminar flow thro- ugh rectangular ducts, International Journal of Thermal Sciences, 8, 16-25 Numerical simulation of the micropolar fluid flow and heat transfer... 569 12. Nakamura S., 1991,Applied Numerical Methods with Software, Prentice-Hall, 442-446 13. Rawool A.S., Mitra S.K., Kandlikar S.G., 2006, Numerical simulation of flow through mi- crochannels with designed roughness,Microfluidics and Nanofluidics, 2, 215-221 14. Roache P.J., Knupp P.M., 1993,CompletedRichardson extrapolation,Communications in Nu- merical Methods in Engineering, 9, 365-374 15. Shangjun Y., Kequn Z., Wang W., 2006, Laminar flow of micropolar fluid in rectangular microchannels,Acta Mechanica Sinica, 22, 403-408 16. SuttonR.S., BartoA.G., 2008,ExactNavier-Stokes solution for pulsatory viscous channel flow with arbitrary pressure gradient, Journal of Propulsion and Power, 24, 1412-1423 Manuscript received June 3, 2013; accepted for print December 6, 2013