Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 52, 4, pp. 1033-1046, Warsaw 2014 NUMERICAL SIMULATION OF UNSTEADY WATER-BASED NANOFLUID FLOW AND HEAT TRANSFER BETWEEN TWO ORTHOGONALLY MOVING POROUS COAXIAL DISKS Kashif Ali, Muhammad Farooq Iqbal, Zubair Akbar, 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 unsteady hydromagnetic (MHD) flow and heat trans- fer characteristics of a viscous incompressible electrically conducting water-based nanofluid (containing Al2O3 nanoparticles) between two orthogonally moving porous coaxial disks with suction. Different from the classical shooting methodology, we employ a combination of a direct and an iterative method (SOR with optimal relaxation parameter) for solving the sparse systems of linear algebraic equations arising from the FD discretization of the linearized self similar nonlinear ODEs. Effects of the governing parameters on the flow and heat transfer are discussed and presented through tables and graphs. The findings of the present investigationmay be beneficial for electronic industry inmaintaining the electronic components under effective and safe operational conditions. Keywords:porousdisks,Al2O3 nanoparticles, viscousdissipation,wall expansion ratio, Joule heating 1. Introduction Theflowbetween two disks havemany important applications in the fields of oceanography, bio- mechanics, rotatingmachinery, heat andmass exchangers, crystal growth processes, lubricants, computer storage devices and viscometries. In thrust bearings, the disks are separated bymeans of a lubricant injected through the disks. Furthermore, fluids with polymer additives have been used as improved lubricating oils in the modern lubrication technology (Connr et al., 1968). In general, it iswell known that conventional heat transfer fluids such asmineral oil, ethylene glycol and water have poor heat transfer properties as compared to those of most solids. To improve heat transfer, very small sized solid particles are added to the base fluid. The resulting new kind of the fluid (named as nanofluid)was firstly introduced byChoi (1995). The nanofluid is used to describe a mixture of a liquid and a solid. Such nanofluids have a very high thermal conductivity which is beyond the explanation of any existing theory. Nanofluids are also more stable and pose no additional problems such as sedimentation, additional pressure drop and erosion, due to the nano sized particles. The comprehensive references on nanofluids can be found in the recent book by Das and Choi (2007) and in the papers by Buongiorno (2006), Kakac and Pramuanjaroenkij (2009). Injection or suction of a fluid through a porous bounding heated or cooled surface is of general interest in many practical applications including film cooling, control of the boundary layer etc. Fluid flow and convective heat transfer in rotor-stator systems (the rotor for a rotating disk and the stator for a stationary disk) are of great importance in turbomachinery and power engineering. Many researchers have worked on the problems related to disks with different wall conditions, for example, Usha and Ravindranthe (2001), Banchok et al. (2011), Elcrat (1976) and Rasmussen (1970). 1034 K. Ali et al. To our best knowledge no researcher has yet considered the behaviour of a nanofluid betwe- en two orthogonally moving porous disks. Therefore, in this paper, we consider the unsteady, laminar, incompressible, and two-dimensional flow of a nanofluid between two orthogonally mo- ving porous coaxial disks, with an external magnetic field acting normally. We use a similarity transformation to reduce the governing partial differential equations to a set of nonlinear co- upled ordinary differential equations in the dimensionless form, which are numerically solved by employing an algorithmbased on the quasi-linearization and finite difference discretization. The ease in obtaining the numerical solution using this techniquemakes it superior than the shooting like approach used in our earlier investigations (Ashraf et al., 2009; Ashraf andWehgal (2012). The effects of the governing parameters on the flow and heat transfer aspects of the problem are discussed. 2. Mathematical formulation Weconsider hydromagnetic unsteady laminar viscous flowandheat transfer of an incompressible nanofluid (containing Al2O3 nanoparticles) between two orthogonally moving porous coaxial diskswith suction in thepresenceof a transversally appliedmagnetic field.The inducedmagnetic field is assumed to be negligible as compared with the imposed field. The magnetic 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 small magnetic Reynolds numbers, the magnetic field tends to relax towards a purely diffusive state. It is also assumed that there is no applied polarization and hence no electric field. The disks have the same permeability, move down or up uniformly at a time-dependent rate a′(t), and are distanced by 2a(t). The geometry of the problem suggests that the cylindrical coordinate systemmay be chosen with the origin at themiddle of the two disks.With u and w being the velocity components in, respectively, the r and z direction, the governing equations for the problem, taking the viscous dissipation and Joule heating effects into account, are ∂u ∂r + u r + ∂w ∂z =0 ∂u ∂t +u ∂u ∂r +w ∂u ∂z =− 1 ρnf ∂p ∂r +υnf (∂2u ∂r2 + 1 r ∂u ∂r − u r2 + ∂2u ∂z2 ) − σeB 2 0 ρnf u ∂w ∂t +u ∂w ∂r +w ∂w ∂z =− 1 ρnf ∂p ∂z +υnf (∂2w ∂r2 + 1 r ∂w ∂r + ∂2w ∂z2 ) ∂T ∂t +u ∂T ∂r +w ∂T ∂z =αnf ∂2T ∂z2 + µnf (ρcp)nf (∂u ∂z )2 + 1 (ρcp)nf σeB 2 0u 2 (2.1) where σe is the electrical conductivity, B0 is the strength of themagnetic field, p is the pressure, and T is the temperature, αnf is the thermal diffusivity, ρnf is thedensity, υnf is thekinematics viscosity of the nanofluid, which are given by υnf = µnf ρnf µnf = µf (1−φ)2.5 ρnf =(1−φ)ρf +φρs αnf = knf (ρcp)nf (ρcp)nf =(1−φ)(ρcp)f +φ(ρcp)s knf kf = ks+2kf −2φ(kf −ks) ks+2kf +φ(kf −ks) (2.2) where ρs and ρf are, respectively, the densities of the solid fractions and of the fluid, (ρcp)nf is the heat capacitance of the nanofluid, knf is the effective thermal conductivity of the nanofluid, approximated by theMaxwell-Garnett model. Numerical simulation of unsteady water-based nanofluid flow... 1035 The boundary conditions are z=−a(t) u=0 w=−Aa′(t) and z= a(t) u=0 w=Aa′(t) (2.3) Here A is the measure of wall permeability, and the dash denotes the derivative w.r.t. time t. After eliminating the pressure term from the governing equations, we introduce the following similarity transformation η= z a u=− rνf a2 Fη(η,t) w= 2νf a F(η,t) θ= T −T2 T1−T2 (2.4) The dimensions of νf are [L 2T−1], those of both u and w are [LT−1], and finally [L] is the dimension of each a and r, which when used in Eq. (2.4) give F = aw/(2νf) and Fη =−a 2u/(rνf) as the two dimensionless velocities in the axial and radial direction, respecti- vely. On the other hand, θ(η) being the ratio of two quantities having the same units, is also dimensionless. The transformation given in Eq. (2.4) leads us to υnf υf Fηηηη +α(3Fηη +ηFηηη)−2FFηηη − a2 υf Fηηt − ρf ρnf MFηη =0 θηη − υf αnf (2F −ηα)θη +[(1−φ) −2.5F2ηη +MF 2 η ]EcPr ( kf knf ) − a2 αnf θt =0 (2.5) with boundary conditions η=−1 F =−R Fη =0 and η=1 F =R Fη =0 (2.6) Here T1 and T2 (with T1 >T2) are fixed temperatures of the lower andupperdisks, respectively, α= aa′(t)/υf is the wall expansion ratio, R =Aaa ′/(2υf) is the permeability Reynolds num- ber, M = σeB 2 0a 2/µf is the magnetic parameter, Pr = (µcp)f/kf is the Prandtl number, and Ec =(rυf) 2/[a4(T1 −T2)(cp)f] is the Eckert number. It is worth to note that continuity Eq. (2.1)1 is identically satisfied, which means that the proposed velocity is compatible with (2.1)1 and, therefore, represents the possible fluidmotion. Finally, we set f =F/R, and consider the case (followingMajdalani et al., 2002) when α is a constant, f = f(η) and θ= θ(η), which leads to θt =0 and fηηt =0. Thus we have υnf υf fηηηη +α(3fηη +ηfηηη)−2Rffηηη − ρf ρnf Mfηη =0 θηη − υf αnf (2Rf −ηα)θη +R 2[(1−φ)−2.5f2ηη +Mf 2 η]EcPr ( kf knf ) =0 (2.7) and η=−1 f =−1 fη =0 and η=1 f =1 fη =0 (2.8) 3. Numerical solution An inspection of Eq. (2.7)1 reveals that it can be readily integrated and becomes υnf υf fηηη +fηη(ηα−2Rf)+fη ( 2α+Rfη − ρf ρnf M ) =β (3.1) where β is the constant of integration to be determined. 1036 K. Ali et al. We consider this third order equation as the following system of coupled 1st and 2nd order ODEs fη = q υnf υf qηη +qη(ηα−2Rf)+q ( 2α+Rq− ρf ρnf M ) =β (3.2) We solve the above equations along with Eqs. (2.7)2, subject to the boundary conditions f(±1)=±1 q(±1)= 0 θ(−1)= 1 θ(1)= 0 (3.3) For the numerical solution of the present problem, we first discretize Eqs. (2.7)2 and (3.2)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.3). On the other hand, Eq. (3.2)1 is integrated numerically by employing the 4th order Runge Kuttamethod, after every SOR iteration. It is important tomention that the boundary condition f(1)= 1 is not needed for the above mentioned computational procedure, and is used in finding the constant of integration β which is the only unknown as compared to the three missing initial conditions in the usual shooting approach. Inour earlierwork,weused trial anderrormethod tofind β but someone-dimensional zero findingalgorithmsmay also be employed for the purpose.However, the sensitivity of βwith respect to the governing parametersmakes it difficult to find it, and somemanual 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 the nonlinear ODEs. 3.1. Quasi-linearization method In quasi-linearization, we construct two sequences of vectors {f(k)} and {θ(k)}, which co- nverge to the numerical solutions to Eqs. (2.7)1 and (2.7)2, respectively. To construct {f (k)}, we linearize Eq. (2.7)1 by retaining only the first order terms, as follows. We set G(f,fη,fηη,fηηη,fηηηη)≡ υnf υf fηηηη +α(3fηη +ηfηηη)−2Rffηηη − ρf ρnf Mfηη G(f(k),f(k)η ,f (k) ηη ,f (k) ηηη,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)ηηη −f (k) ηηη) ∂G ∂f (k) ηηη +(f(k+1)ηηηη −f (k) ηηηη) ∂G ∂f (k) ηηηη =0 υnf υf f(k+1)ηηηη +(ηα−2Rf (k))f(k+1)ηηη + ( 3α− ρf ρnf M ) f(k+1)ηη −2Rf (k) ηηηf (k+1) =−2Rf(k)ηηηf (k) (3.4) Now Eqs. (3.4) gives a system of linear differential equations, with fk being the numerical solution vector of the k-th equation. To solve the linear 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)) (3.5) where n is the number of grid points. Thematrices An×n and Bn×1 are initialized as follows Numerical simulation of unsteady water-based nanofluid flow... 1037 A1,1 =1 B1 =−1 A2,1 =−8d1+2h(αη2 −2Rf2)+2(3α−Md2)h 2 A2,2 =14d1 −4(3α−Md2)h 2−2hR(f4−2f3+2f1−f2)−h(αη2 −2Rf2) A2,3 =−8d1−2h(αη2 −2Rf2)+2(3α−Md2)h 2 A2,4 =2d1+h(αη2 −2Rf2) B2 =−2hRf2(−f2+2f1−2f3+f4) Ai,i−2 =2d1−h(αηi −2Rfi) Ai,i−1 =−8d1+2h(αηi −2Rfi)+2(3α−Md2)h 2 Ai,i =12d1 −4(3α−Md2)h 2 −2hR(fi+2−2fi+1+2fi−1−fi−2) Ai,i+1 =−8d1−2h(αηi −2Rfi)+2(3α−Md2)h 2 Ai,i+2 =2d1+h(αηi −2Rfi) Bi =−2hRfi(fi+2−2fi+1+2fi−1−fi−2)                        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 a rapid convergence, while for a fine grid, smaller valuesmay be chosen for better estimation of ωopt. Wemay improve the order of accuracy of the solution by solving the problem again on the grid with a step h/2 and h/4, and then using theRichardson 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) 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 and the heat transfer rates at the disks which are, respectively, proportional to f ′′(−1) and θ′(−1). Due to symmetry of the problem, the results are given only at the lower disk. The parameters for the present study Numerical simulation of unsteady water-based nanofluid flow... 1039 are the Reynolds number R, the magnetic parameter M, the nanoparticle volume fraction parameter φ, the wall expansion ratio α and the Eckert number Ec. It is worth to note that α < 0 or α > 0 according to the case when the disks are approaching each other or moving away, whereas R < 0 for suction. In the present study, we need to consider the effects of the parameters R,M, φ, α,Ec on the following quantities: • Axial and radial velocity profiles • Temperature distribution • Shear stress and heat transfer rate at the disk. As we are considering the water-based nanofluid containing Al2O3 nanoparticles, we take the following values for the study: ρf = 997.1, ρs = 3970, ks = 40, kf = 0.613, (cp)s = 765, (cp)f =4179 and Pr=6.2. It is to be noted that the case φ=0 corresponds to pure water. Table 1 shows the convergence of our numerical results as the step-size decreases, which gives us confidence on our computational procedure. Itmay be shown, as in our previouswork (Ali et al., 2013), that the shear stress and the heat transfer rate at each disk are proportional to the values of f ′′ and θ′, respectively, at η = −1. Therefore, in the Tables 1-6, we have considered the effects of the governing parameters on f ′′(−1) and θ′(−1). It is clear fromTable 2 that the external magnetic field increases the shear stress as well as the heat transfer rate, whether the disks are approaching each other or moving away. Addition of nanoprticles to water results in remarkably increasing the heat transfer rate while decreasing the shear stress for the two cases of α as predicted by Table 3. On the other hand, suction increases both the shear stress and the heat transfer rate when the disks are moving away, but decreases the shear stress while increasing the heat transfer rate for the other case (clear from Table 4). Table 5 predicts that both the shear stress and heat transfer rate increase when the disks are approaching each other, whereas the opposite trend is noted if the disks aremoving away. Table 6 shows that the viscous dissipation always increases the heat transfer rate at the disks. Table 1.Dimensionless velocity f(η) on three grid sizes and extrapolated values for R=−10, φ=0.3, α=−5 and M =1 f(η) η 1st grid 2nd grid 3rd grid Extrapolated (h=0.02) (h=0.01) (h=0.005) values 0 0 0 0 0 0.2 0.3054554 0.3054592 0.3054602 0.3054605 0.4 0.5806770 0.5806883 0.5806911 0.5806920 0.6 0.8011234 0.8011416 0.8011461 0.8011476 0.8 0.9470534 0.9470702 0.9470744 0.9470758 Table 2. Effect of the magnetic parameter M on f ′′(−1) and θ′(−1) for φ = 0.2, R = −10 and Ec=0.1 M α=2 α=−2 f ′′(−1) θ′(−1) f ′′(−1) θ′(−1) 0 1.8580 2.9180 2.3175 5.4619 2 1.9009 3.0519 2.3685 5.6988 4 1.9439 3.1893 2.4194 5.9406 6 1.9871 3.3300 2.4704 6.1873 8 2.0303 3.4740 2.5213 6.4388 1040 K. Ali et al. Table 3. Effect of the nanoparticle volume fraction parameter φ on f ′′(−1) and θ′(−1) for M =2,R=−10 and Ec=0.1 φ α=2 α=−2 f ′′(−1) θ′(−1) f ′′(−1) θ′(−1) 0 1.9200 1.6794 2.3936 3.1664 0.05 1.9103 1.9174 2.3823 3.6112 0.1 1.9044 2.2135 2.3749 4.1606 0.15 1.9013 2.5839 2.3704 4.8430 0.2 1.9009 3.0519 2.3685 5.6988 Table 4.Effect of thepermeabilityReynoldsnumber R on f ′′(−1) and θ′(−1) for and Ec=0.1 R α=2 α=−2 f ′′(−1) θ′(−1) f ′′(−1) θ′(−1) −2 1.7950 0.4395 3.5007 3.0065 −4 1.8234 1.0183 2.9193 3.8376 −6 1.8573 1.6661 2.6265 4.4320 −8 1.8826 2.3498 2.4660 5.0472 −10 1.9009 3.0519 2.3685 5.6988 Table 5. Effect of the wall expansion ratio α on f ′′(−1) and θ′(−1) for φ = 0.2, R = −10, M =2 and Ec=0.1 α f ′′(−1) θ′(−1) −5 2.8351 9.5692 −2 2.3685 5.6988 −1 2.2370 4.8414 1 2.0041 3.5434 2 1.9009 3.0519 5 1.6355 1.9997 Table 6.Effect of the Eckert number Ec on θ′(−1) for R=−10, φ=0.2 and M =2 Ec θ′(−1) α=2 α=−2 0 0 0 0.1 3.0519 5.6988 0.2 6.1039 11.3975 0.3 9.1558 17.0963 0.4 12.2077 22.7951 The streamlines for the present problem are shown in Fig. 1b. Figs. 2 and 3 show that, for α> 0 or α< 0, the externalmagnetic field tends to decrease both the axial and radial velocities in the region in the middle of the two disks. Thus, for this region, the magnetic field acts like a drag force called the Lorentz force which decreases the fluid velocity. This results in generation of heat energy which remarkably increases the fluid temperature as shown in Figs. 2c and 3c. The influence of the nanoparticle volume fraction parameter on the velocity and temperature profiles is similar to that of the external magnetic field (Figs. 4 and 5) for the two cases of α. The effect of R on the velocity and temperature distribution is the same as that of M for α > 0, whereas the opposite trend may be observed for α < 0 (Figs. 6 and 7). As the values Numerical simulation of unsteady water-based nanofluid flow... 1041 of α are varied from negative to positive, the axial velocity component increases, whereas the radial velocity increases in the middle of the region between the disks and decreases near the disks (Figs. 8a and 8b).When the disks are approaching each other, an increase in the expansion ratio α results in increasing the temperature distribution across the disks. But in the other case, the temperature profiles rise only in the middle of the disks (Fig. 8c). Finally, Figs. 9a and 9b show that the viscous dissipation remarkably increases the temperature distribution across the disks, nomatter if the disks are approaching or receding. Fig. 1. (a) Physical model of the problem; (b) streamlines for the problem for R=−10, φ=0.3, α=−5 and M =1 Fig. 2. (a) Axial velocity, (b) radial velocity and (c) temperature profiles for R=−2, φ=0.1, α=2 and Ec=0.1 1042 K. Ali et al. Fig. 3. (a) Axial velocity, (b) radial velocity and (c) temperature profiles for R=−2, φ=0.1, α=−2 and Ec=0.1 Fig. 4. (a) Axial velocity, (b) radial velocity and (c) temperature profiles for M =1, α=1,R=−10 and Ec=0.1 Numerical simulation of unsteady water-based nanofluid flow... 1043 Fig. 5. (a) Axial velocity, (b) radial velocity and (c) temperature profiles for M =1, α=−1,R=−10 and Ec=0.1 Fig. 6. (a) Axial velocity, (b) radial velocity and (c) temperature profiles for M =1, φ=0.1, α=5 and Ec=0.1 1044 K. Ali et al. Fig. 7. (a) Axial velocity, (b) radial velocity and (c) temperature profiles for M =1, φ=0.1, α=−5 and Ec=0.1 Fig. 8. (a) Axial velocity, (b) radial velocity and (c) temperature profiles for M =1, φ=0.2,R=−3 and Ec=0.1 Numerical simulation of unsteady water-based nanofluid flow... 1045 Fig. 9. Temperature profiles for: (a) M =1, α=5,R=−10 and φ=0.2, (b) M =1, α=−5, R=−10 and φ=0.2 5. Conclusion In this paper, we have numerically studied how the governing parameters, namely themagnetic parameter M, theReynolds number R, the nanoparticle volume fraction parameter φ, the wall expansion ratio α and the Eckert number Ec affect the flow and heat transfer characteristics of the unsteady, laminar, incompressible and two-dimensional flow of a water-based nanofluid (containing the Al2O3 nanoparticles) between two orthogonally moving porous coaxial disks with suction. When the disks are moving away (α> 0): • shear stress at the disks increases with M and R, whereas the opposite effect is observed for φ and α, • heat transfer rate at the disks increases with M,R and φ. For the approaching disks (α< 0): • shear stress at the disks increases with M and α, whereas the opposite effect is observed for φ and R, • heat transfer rate at the disks increases with M,R, α and φ. All the governing parameters significantly increase the temperature distribution in the middle of the region between the approaching or receding disks. Finally, it is important tomention that the heat transfer characteristics of the nanofluids are not yet fully explored, and that is why complementary works are necessary to identify new and unique applications of these fluids. Acknowledgements The authors are extremely grateful to theHigher EducationCommission of Pakistan for the financial support to carry out this research. The authors are also extremely grateful to the learned reviewer for useful comments improving quality of this paper. References 1. Ali K., Ashraf M., Jameel N., 2013, Numerical simulation of MHDmicropolar fluid flow and heat transfer in achannelwith shrinkingwalls,Canadian Journal ofPhysics, 10.1139/cjp-2013-0324 2. AshrafM.,KamalM.A., SyedK.S., 2009,Numerical simulation of amicropolar fluid between a porousdisk andanon-porousdisk,Journal of AppliedMathematics andModelling, 33, 1933-1943 3. Ashraf M., Wehgal A.R., 2012,MHD flow and heat transfer of micropolar fluid between two porous disks,Applied Mathematics and Mechanics (English Edition), 33, 51-64 1046 K. Ali et al. 4. Banchok N., Ishak A., Pop I., 2011, Flow and heat transfer over a rotating porous disk in a nanofluid,Physica B, 406, 1767-1772 5. Buongiorno J., 2006, Convective transport in nanofluids,ASME Journal of Heat Transfer, 128, 240-250 6. Choi S.U.S., 1995, Enhancing thermal conductivity of fluids with nano particles, [In:] Develop- ments and Applications of Non-Newtonian Flows, Siginer D.A.,Wang H.P. (Eds.), 31, 99-105 7. Connor J.J., Boyd J., Avallone E.A., 1968, Standard Handbook of Lubrication Engineering, McGraw-Hill, NewYork 8. Das S.K., Choi S.U.S., 2007,Nanofluids: Science and Technology, Wiley, New Jersey 9. Deuflhard P., 1983,Order and step-size control in extrapolationmethods,Numerical Mathema- tics, 41, 399-422 10. Elcrat A. R., 1976, On the radial flow of a viscous fluid between porous disks, Archive for Rational Mechanics and Analysis, 61, 91-96 11. Kakac S., Pramuanjaroenkij A., 2009, Review of convective heat transfer enhancement with nanofluids, International Journal of Heat and Mass Transfer, 52, 3187-3196 12. Majdalani J., Zhou C., Dawson C.A., 2002, Two dimensional viscous flows between slowly expanding or contracting walls with weak permeability, Journal of Biomechanics, 35, 1399-1403 13. Nakamura S., 1991,Applied Numerical Methods with Software, Prentice-Hall, 442-446 14. Rasmussen H., 1970, Steady viscous flow between two porous disks, Zeitschrift für Angewandte Mathematik und Physik, 21, 187-195 15. Roache P.J., Knupp P.M., 1993,CompletedRichardson extrapolation,Communications in Nu- merical Methods in Engineering, 9, 365-374 16. UshaR.,RavindrantheR., 2001,Numerical studyof filmcooling on rotatingdisk, International Journal of Nonlinear Mechanics, 36, 147-154 Manuscript received January 13, 2014; accepted for print June 24, 2014