SQU Journal for Science, 2016, 21(1), 48-63 © 2016 Sultan Qaboos University 48 Effects of Second-Order Slip and Viscous Dissipation on the Analysis of the Boundary Layer Flow and Heat Transfer Characteristics of a Casson Fluid Mohammad M. Rahman 1* and Ioan Pop 2 1 Department of Mathematics and Statistics, College of Science, Sultan Qaboos University, P.O. Box: 36, PC 123, Al-Khod, Muscat, Sultanate of Oman. 2 Department of Mathematics, Faculty of Mathematics and Computer Science, Babeş-Bolyai University, 400084 Cluj- Napoca, Romania. *Email: mansur@squ.edu.om. ABSTRACT: The aim of the present study is to analyze numerically the steady boundary layer flow and heat transfer characteristics of Casson fluid with variable temperature and viscous dissipation past a permeable shrinking sheet with second order slip velocity. Using appropriate similarity transformations, the basic nonlinear partial differential equations have been transformed into ordinary differential equations. These equations have been solved numerically for different values of the governing parameters namely: shrinking parameter , suction parameter ,s Casson parameter , first order slip parameter ,a second order slip parameter ,b Prandtl number Pr, and the Eckert number Ec using the bvp4c function from MATLAB. A stability analysis has also been performed. Numerical results have been obtained for the reduced skin-friction, heat transfer and the velocity and temperature profiles. The results indicate that dual solutions exist for the shrinking  0  surface for certain values of the parameter space. The stability analysis indicates that the lower solution branch is unstable, while the upper solution branch is stable and physically realizable. In addition, it is shown that for a viscous fluid     a very good agreement exists between the present numerical results and those reported in the open literature. The present results are original and new for the boundary-layer flow and heat transfer past a shrinking sheet in a Casson fluid. Therefore, this study has importance for researchers working in the area of non-Newtonian fluids, in order for them to become familiar with the flow behavior and properties of such fluids. Keywords: Boundary layer; Casson fluid; Dual solutions; Shrinking surface; Second-order slip; Stability analysis. تأثير اإلوزالق مه انذرجة انثاوية و انتبذد انهزج عهى تحهيم انتذفك عبر انطبقة انمتاخمة و خصائص اإلوتقال انحراري نسائم كاسون محمذ مىصور رحمه و آيون بوب ٌٙذف ٘زا اٌبحث إٌى دساست اٌتذفك اٌّطشد عبش اٌطبمت اٌّتاخّت ٚخصائض إٔتماي اٌحشاسي ٌسائً واسْٛ ِع ٚخٛد دسخاث حشاسة :انمهخص ٔزالق ِٓ اٌذسخت اٌثأٍت. ٚتعتّذ طشٌمت اٌحً عٍى استخذاَ تحٌٛالث اٌتشابٗ إِتغٍشة ٚتبذد ٌزج حٛي سطح رٚ ٔفارٌت ٚلابً ٌإلٔىّاش ٚفً ظً ٚخٛد ٚاٌتً ٌتُ ِٓ خالٌٙا تحًٌٛ اٌّعادالث األساسٍت فً إٌّٛرج اٌشٌاضً اٌّستخذَ ٚاٌتً ً٘ عباسة عٓ ِعادالث تفاضٍٍت خزئٍت غٍش خطٍت إٌّاسبت ٌداد حٍٛي عذدٌت ٌالحتىان اٌسطحً إ. ٚتُ "MATLABإٌى ِعادالث تفاضٍٍت اعتٍادٌت ٚاٌتً ٌتُ حٍٙا بطشق اٌتحًٍٍ اٌعذدي باستخذاَ بشٔاِح "  ماي اٌحشاسي ٚسشعت تذفك ٚدسخاث حشاسة اٌسائً باستخذاَ لٍُ ِختٍفت ٌعذة عٛاًِ فٍزٌائٍت ًٚ٘ ِعاًِ االٔىّاشٚاالٔت ٚ ِعاًِ االِتصاص أٚ ، وّا تُ اٌمٍاَ بتحًٍٍ ثباث Ecٚ إوٍشث Prٚأعذاد باسٔذتً bٚاٌثأٍت aٚ ِعاِالث اإلٔزالق ِٓ اٌذسخت األٌٚى ِٚعاًِ واسْٛ sاٌشفظ ( ٚرٌه ٌمٍُ ِحذد ٌٍّعاِالث اٌفٍزٌائٍت، وّا 0اٌحٍٛي اٌعذدٌت أٌضا. ٚلذ أشاسث إٌتائح إٌى ٚخٛد حٍٛي ِزدٚخت فً حاي األسطح إٌّىّشت ) ت إٌى ٚخٛد أظٙشث ٔتائح تحًٍٍ االستمشاس إٌى عذَ استمشاس اٌدزء اٌسفًٍ ِٓ اٌحً بٍّٕا اٌدزء اٌعٍٛي ِٓ اٌحً ِستمش ٌّٚىٓ تحمٍمٗ فٍزٌائٍا، إضاف )تعٍك باٌسائً اٌٍزج تٛافك بٍٓ ٔتائح اٌذساست اٌحاٌٍت ٚاٌذساساث اٌسابمت فٍّا ٌ    . إْ إٌتائح اٌحاٌٍت تعذ لٍّت ٚخذٌذة فً دساست اٌتذفك ( ِداي اٌّطشد عبش اٌطبمت اٌّتاخّت ٚخصائض االٔتماي اٌحشاسي ٌسائً واسْٛ حٛي سطح لابً ٌالٔىّاش. ٚتعذ ٘زٖ اٌذساست راث أٍّ٘ت ٌٍباحثٍٓ فً ت سٍٛن ٚخٛاص ٘زٖ إٌٛعٍت ِٓ اٌسٛائً.اٌسٛائً غٍش ٍٔٛتٍٍٕت ٌّعشف سائً واسْٛ، اٌطبمت اٌّتاخّت، إٔىّاش األسطح، اإلٔزالق ِٓ اٌذسخت اٌثأٍت، حٍٛي ِزدٚخت، تحًٍٍ االستمشاس.: انكهمات انمفتاحية EFFECTS OF SECOND-ORDER SLIP AND VISCOUS DISSIPATION 49 1. Introduction n the past several decades, there has been an increasing interest in the flows of Newtonian and non-Newtonian fluids over stretching/shrinking sheets because of their applications in processing industries such as paper production, hot rolling, wire drawing, glass-fiber production, aerodynamic extrusion of polymer fiber extruded continuously from a dye with a tacit assumption that the fiber is inextensible, etc. The cooling of a long metallic wire in a bath (an electrolyte) is another physical situation belonging to this category. Glass blowing, continuous casting, and spinning of fibers also involve the flow due to a stretching surface. During its manufacturing process a stretched sheet interacts with the ambient fluid thermally and mechanically. The thermal interaction is governed by the surface heat flux. This surface heat flux can either be prescribed or be the output of a process in which the surface temperature distribution has been prescribed. However, in real life situations one encounters a boundary layer flow over a non-linear stretching surface. For example, in a melt-spinning process, the extrudate is stretched into a filament while it is drawn from the dye. Finally, this surface solidifies while it passes through an effectively controlled cooling system in order for the final product to achieve top quality. Other examples include drawing of copper wires, continuous stretching of plastic films and artificial fibers, hot rolling, glass-fiber, metal extrusion, metal spinning, etc. (see Sparrow and Abraham [1], and Abraham and Sparrow [2]). The pioneering work on the steady boundary layer flow due to a linear stretching sheet was done by Crane [3]. Thereafter, various aspects of this problem have been studied extensively in Newtonian and non-Newtonian fluids. However, the extension of the theory of Newtonian fluids to the theory of non-Newtonian fluids has been proved to be not so straightforward (see, for example, Skelland [4], Denn [5], Rajagopal et al. [6], Bird et al. [7] and Slattery [8]). It is now generally recognized that, in real industrial applications, non-Newtonian fluids are more appropriate than Newtonian fluids. These fluids have wide-ranging industrial applications, such as in the design of thrust bearings and radial diffusers, drag reduction, transpiration cooling, thermal oil recovery, etc. In certain polymer processing applications, one deals with the flow of a second-order (viscoelastic) fluid over a stretching surface. Such fluids are referred to as fluids of the differential type, that is, fluids whose stress is determined by the Rivlin–Ericksen tensors (see Rivlin and Ericksen [9]), or fluids of the rate type, such as the Oldroyd-B fluid (see Oldroyd [10]), 3 rd grade fluid (see Hayat et al. [11]) etc. Polymers mixed in Newtonian solvents and polymer melts, such as high-viscosity silicone oils or molten plastics, are examples of such fluids. But in practice, many materials, like melts, muds, printing ink, condensed milk, glues, soaps, shampoos, sugar solution, paints, tomato paste, etc., show various characteristics which are not properly understandable using Newtonian theory. Hence, to analyse such fluids it is essential to utilise non-Newtonian fluid mechanics. However, the main difficulty is to make a single constitutive equation which covers all the properties of such non-Newtonian fluids. Because of this, numerous non-Newtonian fluid models have been proposed in the literature. In nature, some non-Newtonian fluids behave like elastic solids, i.e. with small shear stress no flow occurs. Casson fluid is one such fluid. So, for flow to occur, the shear stress magnitude of a Casson fluid needs to exceed the yield shear stress, as otherwise the fluid behaves as a rigid body. This type of fluid can be categorised as a purely viscous fluid with high viscosity. Several models have been suggested for non-Newtonian fluids, with their constitutive equations varying greatly in complexity. Some authors (Fredrickson [12]) studied Casson fluid for the flow between two rotating cylinders. Eldabe and Elmohands [13] investigated the steady-flow behavior of a viscoelastic flow through a channel bounded by two permeable parallel plates, and Boyd et al. [14] discussed steady and oscillatory blood flow taking into account Casson fluid. Mernone et al. [15] described the peristaltic flow of a Casson fluid in a two- dimensional channel. Mustafa et al. [16] studied the unsteady boundary layer flow and heat transfer of a Casson fluid over a moving flat plate with a parallel free stream and solved this problem analytically using the homotopy analysis method. Pramanik [17] studied Casson fluid flow and heat transfer past an exponentially porous stretching surface in the presence of thermal radiation. The objective of the present study is to analyze the boundary layer development (both momentum and thermal energy) of a Casson fluid past a permeable shrinking sheet with second-order slip velocity and viscous dissipation. Applying similarity transformations, the basic nonlinear partial differential equations are transformed into ordinary ones which are then solved numerically for different values of the governing parameters. A stability analysis is carried out for the acceptance of the physically realizable solutions. Numerical results are obtained for the reduced skin-friction coefficient, rate of heat transfer and the velocity and temperature distributations. I MOHAMMAD M. RAHMAN AND IOAN POP 50 Figure 1. The geometry of the problem of a shrinking surface. 2. Basic Equations Consider the steady two-dimensional flow of a Casson fluid past a permeable shrinking sheet with the velocity )(xu w , while the second-order wall velocity slip is )( slip xu . The sheet is situated at 0y , the flow being confined to 0y , where y is the coordinate measured in the normal direction to the surface and the coordinate x is measured along the surface of the sheet as shown in Figure 1. It is assumed that the temperature of the surface of the sheet is ( ), w T x while the constant temperature of the ambient (inviscid) fluid is ,T  where   TxT w )( (heated sheet). We assume also that the rheological equation of state for an isotropic and incompressible flow of a Casson fluid can be written as, nnn /1/1 0 /1   (1) or, (see Nakamura and Sawada [18]) ji n n y Bji e p 2 2 /1                    (2) where  is the dynamic viscosity, B  is plastic dynamic viscosity of the non- Newtonian fluid, y p is yield stress of fluid,  is the product of the component of deformation rate with itself, namely, ijij ee , and ij e is the th ji ),( component of the deformation rate, and that the constant n takes the value 1n in many practical applications. Under these conditions, and along with the assumption that the viscous dissipation term in the energy equation is taken into consideration, the boundary layer equations which govern this problem are 0 u v x y       (3) 2 2 (1 1 / ) u u u u v x y y            (4) 2 2 2 (1 1 / ) p T T T u u v x y y C y                      (5) subject to the initial and boundary conditions 0 slip , ( ) ( ), ( ) at 0 0, as w w v v u u x u x T T x y u T T            (6) u w T w T y x  0 Boundary layer v 0 u w EFFECTS OF SECOND-ORDER SLIP AND VISCOUS DISSIPATION 51 where u and v are the velocity components along x and y axes, T is the fluid temperature,  is the thermal diffusivity of the fluid,  is the kinematic viscosity, p C is the specific heat at constant pressure,  is the constant stretching ( 0 ) or shrinking ( 0 ) parameter, 2 / B c y p   is the non-Newtonian (Casson) parameter, and where we assume ( ) w u x cx and 2 0 ( ) w T x T T x    , )0(0 T to be the temperature characteristic of the sheet. Further, following Wu [19] (used also by Fang et al. [20], Mahmood et al. [21]), the wall second-order slip velocity )( slip xu , is given by 2 2 2 2 2 2 3 3 1 1 24 2 ( ) (1 ) = Aslip 2 2 23 2 4 l l u u u u u x l l B K y yKn y yn                              (7) Here )1,/1min( n Kl  and  is the momentum accommodation coefficient with 10  . Based on the definition of l , it is seen that for any given value of n K , we have 10  l . Since the molecular mean free path  is always positive it results in B being a negative number. It is important to mention that Wu’s [19] second order slip velocity model is valid for any arbitrary Knudsen number and it matches better with the Fukui–Kaneko results which are based on the direct numerical simulation of the linearized Boltzmann equation [22]. We look for a similarity solution of Equations (4)-(7) of the form: ycTTTTfxc w  /),/()()(),(   (8) where  is the stream function, which is defined in the classical form as yu  / and xv  / . Thus, we have )(),('  fcvfxcu  (9) where prime denotes differentiation with respect to  . Substituting (7) into Eqs. (4) and (5), we obtain the following ordinary (similarity) equations 2 (1 1 / ) ''' '' ' 0f f f f    (10) 21 '' ' 2 ' (1 1 / ) '' 0 Pr f f Ec f        (11) subject to the boundary conditions     as0)(,0)(' 1)0(),0(''')0('')0(',)0( f fbfafsf (12) where cvs / 0  is the constant parameter of suction ( 0s ) or injection ( 0s ), a is the first order velocity slip parameter with 0/  cAa , b is the second order slip velocity with 0/  cBb ,  /Pr  is the Prandtl number and ( ) / ( ( ) )w p wEc u x C T x T    is the constant Eckert number. It is worth pointing out that for  , Eq. (10) reduces to Eq. (7) for a viscous (Newtonian) fluid, as shown by Fang et al. [20]. The physical quantities of interest are the skin friction coefficient f C and the local Nusselt number xNu , which can be easily shown to be given by 1/ 2 1/ 2 Re (1 1 / ) ''(0), Re '(0) x f x x C f Nu       (13) where Re ( ) / x w u x x  is the local Reynolds number. 3. Stability Analysis Solving Eqs. (10) and (11) with the boundary conditions (12), it has been shown that dual (upper and lower branch) solutions exist for the case of a shrinking sheet )0(  . Thus, in order to see which solution is stable, and that it is physically realizable in practice, a stability analysis of these solutions is necessary. In this respect, we write Eqs. (1)-(3) corresponding to the unsteady flow and heat transfer as 0 u v x y       (14) MOHAMMAD M. RAHMAN AND IOAN POP 52 2 2 (1 1 / ) u u u u u v t x y y               (15) 2 2 2 (1 1 / ) p T T T T u u v t x y y C y                         (16) subject to the initial and boundary conditions 0 0 slip 0 : , 0, for any , 0 : , ( ) ( ), ( ) at 0 0, as w w t v v u T T x y t v v u u x u x T T x y u T T                  (17) where t denotes time. Following Weidman et al. [23], Rosca and Pop [24-25], and Rahman et al. [26-28], a dimensionless time variable  has to be introduced. The use of  is associated with an initial value problem and this is consistent with the question of which solution will be obtained in practice (physically realizable). Thus, the new variables for the unsteady problem are ( , ), ( , ) ( ) / ( ), / , w c x f T T T T c y c t                  (18) Thus, Eqs. (14)-(16) can be written as 2 3 2 2 3 2 (1 1 / ) 0 f f f f f                          (19) 2 2 2 2 2 1 2 (1 1 / ) 0 Pr f f f Ec                                 (20) subject to the boundary conditions     as0),(,0),(' 1),0(),,0('''),0(''),0(',),0( f fbfafsf (21) To test the stability of the steady flow solution )()( 0  ff  and )()( 0   satisfying the boundary-value problem (19)-(21), we write (see Weidman et al. [23] or Rahman et al. [26-28]) ),()(),(),,()(),( 00   GeFeff   (22) where  is an unknown eigenvalue parameter, and ),( F and ),( G are small relative to )(0 f and )( 0  Introducing (22) into Equations. (19) and (20), we get the following linearized problem 3 2 2 0 0 03 2 (1 1 / ) ( 2 ') '' 0 F F F f f f f F                       (23) 2 0 0 02 2 0 2 1 ( 2 ') 2 Pr 2 (1 1 / ) '' 0 G G F f f G F Ec f                            (24) along with the boundary conditions 2 3 2 3 (0, ) 0, (0, ) (0, ) (0, ), (0, ) 0 ( , ) 0, ( , ) 0 as F F F F a b G F G                               (25) The solution )()( 0  ff  and )()( 0   of the steady equations (10) and (11) is obtained by setting 0 . Hence )( 0 FF  and )( 0 GG  in (23) and (24) identify initial growth or decay of the solution (22). With respect to this, we have to solve the linear eigenvalue problem 0''')'2(''''')/11( 0000000  FfFfFfF  (26) EFFECTS OF SECOND-ORDER SLIP AND VISCOUS DISSIPATION 53 0 0 0 0 0 0 0 0 0 1 '' ' ( 2 ') 2 ' 2 (1 1 / ) '' '' 0 Pr G f G f G F Ec f F         (27) along with the boundary conditions 0 0 0 0 0 0 0 (0) 0, '(0) ''(0) '''(0), (0) 0 '( ) 0, ( ) 0 as F F a F b F G F G           (28) It should be stated that for particular values of  ,  , a , b , Ec , and Pr the stability of the corresponding steady flow solution )( 0 f and )( 0  is determined by the smallest eigenvalue  . As has been suggested by Harris et al. [29], the range of possible eigenvalues can be determined by relaxing a boundary condition on )( 0 F or )( 0 G . For the present problem, we relax the condition that 0 ( ) 0G   as  and for a fixed value of  we solve the system (26)-(27) along with the new boundary condition 0 (0) 1G  . 4. Numerical Technique Following Rahman et al. [26-28], the ordinary differential equations (10)-(11) subject to the boundary conditions (12) are solved numerically using the function bvp4c from the very robust computer algebra software, MATLAB. The function bvp4c requires writing equations (10)-(11) as first order differential equations by introducing new variables: one for each variable in the original problem plus one for each derivative up to the highest order derivative minus one. It then implements a collocation method for the solution of the following boundary value problem ( , )x w t x  , a t b  (29) subject to the two-point boundary conditions ( ( ), ( )) 0bc x a x b  . (30) The approximate solution ( )t is a continuous function that is a cubic polynomial on each subinterval 1 [ , ] n n t t  of the mesh 0 1 2 .......... N a t t t t b     . It satisfies the boundary conditions ( ( ), ( )) 0bc a b   (31) and it also satisfies the following differential equations (collocates) at both ends and mid-point of each subinterval ( ) ( , ( )) n n n t w t t   (32) 1 1 1 (( ) / 2) (( ) / 2, (( ) / 2)) n n n n n n t t w t t t t          (33) 1 1 1 ( ) ( , ( )) n n n t w t t       (34) These conditions result in a system of nonlinear algebraic equations for the coefficients defining ( )t , which are solved iteratively by linearization. Here ( )t is a fourth order approximation to an isolated solution ( )x t , i.e. 4 || ( ) ( ) ||x t t K h  , where K is the maximum of the step sizes 1n n n h t t    and K is a constant. For such an approximation, the residual ( )r t is the ordinary differential equation and is defined by ( ) ( ) ( , ( ))r t t w t t   . (35) In this approach mesh selection and error control are based on the residual of the continuous solution. The relative error tolerance was set to 7 10  . The condition   needs to be replaced by a suitable finite value of  , say   . We started the computation at a small value, for example, 5 , then subsequently increased the value of  until the boundary conditions were verified. In using this method, we chose a suitable finite value of  , namely 20   for the upper branch (first) solution and   in the range 40-60 for the lower branch (second) solution. The present problem may have more than one solution. Thus the bvp4c function requires an initial guess of the desired solution for the system (10)-(11). The guess should satisfy the boundary conditions and reveal the characteristics of the solution. The bvp4c method always converges to the first solution even for poor guesses of the initial conditions. Thus, determining an initial guess for the first (upper branch) solution is not difficult. On the other hand, it is very difficult to come up with a sufficiently good guess for the second (lower branch) solution of the system (10)-(11). To overcome this difficulty, we need to start with a set of parameter values for which the problem is easy to solve. Then, we use the obtained result as the initial guess for the solution of the problem with small variation of the parameters. This is repeated until we reach the right values of the parameters. This technique MOHAMMAD M. RAHMAN AND IOAN POP 54 is called continuation. Examples of solving boundary value problems by bvp4c can be found in the book by Shampine et al. [30] or through an online tutorial by Shampine et al. [31]. The numerical simulations are carried out for various values of the physical parameters such as Casson parameter (  ), suction parameter ( s ), shrinking parameter (  ), first order slip parameter ( a ) and second order slip parameter ( b ), Eckert number ( Ec ) and Prandtl number ( Pr ). Because of the almost complete lack of experimental data, the choice of the values of the parameters is dictated by the values chosen by previous investigators. The value of the Prandtl number is set equal to 1 throughout the paper unless otherwise specified. The values of the other parameters are mentioned in the description of the respective figures. It is worth mentioning that for a viscous fluid (   ) and for 1  (stretching sheet), in the absence of fluid suction 0s  and second order slip parameter 0b  , equation (10) 2 0f ff f     subject to the boundary conditions (0) 0f  , (0) 1 (0)f af   , ( ) 0f    matches with that of Wang [32-33], Anderson [34], and Sahoo and Do [35]. Wang [32-33] and Anderson [34] obtained an exact solution of (9) subject to the above-mentioned boundary conditions. To test the accuracy of the current numerical solution, the values of (0)f  and ( )f  are compared with analytical values reported by Wang [32-33] and Sahoo and Do [35] in Table 1. This table shows that the numerical values produced by the current code and the exact analytical solutions reported by Wang [32-33] and Shahoo and Do [35] are in very good agreement. This gives us confidence to use the present code. It can further be mentioned that for a viscous fluid (    ) Eq. (10) 2 0f f f f     along with the boundary conditions (0)f s , (0) (0) (0)f af bf     , ( ) 0f   as   exactly matches with the corresponding Eq. (10) and boundary conditions (12) of Rosca and Pop [24] when the buoyancy force in their model is neglected. We also mentioned that for 1   (shrinking sheet), the above-stated equation together with the corresponding boundary conditions also matches with Eqs. (7)-(8) of Fang et al. [20]. Table 2 shows the comparison of the values of (0)f  for 1a  and several values of s and b when 1   with those reported by Fang et al. [20], and Rosca and Pop [24]. In fact, the results show excellent agreement among the data, thus giving us confidence to use the present MATLAB code. 5. Results and Discussion The numerical simulation of Eqs. (10) to (11) subject to the boundary conditions (12) are carried out for various values of the physical parameters , ,s , ,a ,b Ec and Pr for obtaining the condition under which the dual (upper and lower branch) solutions for the steady flow of a Casson fluid over a shrinking surface may exist. Miklavčič and Wang [36] have studied the steady viscous (Newtonian) fluid flows over a permeable linearly shrinking surface and have shown that suction at the wall will generate dual solutions only when the suction parameter s is greater than or equal to 2. Table 1. Comparison of results (0)f  and ( )f  with the first order slip parameter a when   , 0,s  and 1  . a (0)f  ( )f  Current result Sahoo and Do [35] Wang [32] Current result Sahoo and Do [35] Wang [32] Wang [33] 0.0 0.2 0.3 0.5 1.0 2.0 3.0 5.0 10 20 1.000000 0.776377 0.701548 0.591195 0.430159 0.283979 0.214054 0.144840 0.081242 0.043788 1.001154 0.774933 0.699738 0.589195 0.428450 0.282893 0.213314 0.144430 0.081091 0.043748 1.0 - 0.701 - 0.430 0.284 - 0.145 - 0.0438 0.999973 0.919076 0.888544 0.839278 0.754859 0.657267 0.598141 0.524529 0.431842 0.349999 1.001483 0.919010 0.888004 0.838008 0.752226 0.652253 0.590892 0.513769 0.413655 0.322559 1.0 - 0.887 - 0.748 0.652 - 0.514 - 0.332 1.0 - - 0.8393 0.7549 - 0.5982 - 0.4331 - EFFECTS OF SECOND-ORDER SLIP AND VISCOUS DISSIPATION 55 Table 2. Comparison of (0)f  for 1a  , 2, 3s  and 1, 2b    when 1   . s a b Present study Rosca and Pop [24] Fang et al. [20] Lower branch Upper branch Lower branch Upper branch Lower branch Upper branch 2 2 3 3 1 1 1 1 -1 -2 -1 -2 0.25659275 0.22573384 0.20223817 0.18809181 0.29054789 0.18465688 0.23201653 0.1369054 0.2565 0.2257 0.2022 0.1868 0.2905 0.1846 0.2320 0.1369 0.2565 0.2257 0.2022 0.1868 0.2905 0.1847 0.2317 0.1371 Figure 2. Variations of (0)f  for different values of  at 2  when 1,a  1b   , 2.5,s  1,Ec  and Pr 1 . Figure 3. Variations of (0) for different values of  at 2  when 1a  , 1b   , 2.5s  , 1Ec  , and Pr 1 . In Figures 2 to 7, we have investigated the variation of the reduced skin friction coefficient   1/ 2 / ( 1) Re ''(0) x f C f    and the reduced local Nusselt number (or heat transfer from the surface of the sheet) 1/ 2 Re '(0) x x Nu    of a Casson fluid for different values of  , s and  respectively, keeping the values of the other parameters 1a  , 1b   , 1Ec  , and Pr 1 fixed. From these figures we see that the number of solutions depends on the shrinking parameter  , suction parameter s , and Casson parameter  . In Figures 2 and 3 we have identified two critical  ’s say s  and c  such that 0 c s    where s  is the critical value of  for the lower branch and c  is the critical value of  for the upper branch in which the solution exists. That is, for a Casson fluid flow over a shrinking sheet, there exist dual solutions when 0 s    , one solution which is the upper branch when , c s     and no solution when c  . Thus, for c  the full Navier-Stokes equations and energy equation need to be solved. Figure 2 also shows that values of (0)f  for both the upper and lower branches decrease with the increase of  when other parameter values are fixed. For the upper branch solution these values are higher than the corresponding lower branch solution. Now the question is which of these solutions is physically acceptable. From the stability analysis it is found that the upper branch solution is stable and physically acceptable whereas the lower solution branch is unstable and physically unacceptable. Table 3 presents the smallest eigenvalues  for the upper branch solution at several values of ,a b and . This table shows that smallest eigenvalues for the upper branch solution are positive, hence the perturbed part F of the solution f will be diminished and converged to the steady state solution 0 f when   as can be seen from Eq. (22). Thus, the upper branch solution is stable and physically acceptable. On the other hand the smallest eigenvalues for the lower branch solution are negative which indicates that the perturbed part F of the -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 -0.2 0 0.2 0.4 0.6 0.8 1 1.2  f   ( 0 ) [ s , f  (0)] = [-1.5357, 0.440185] [ c , f  (0)] = [-3.0993, 1.028130] upper branch solution lower branch solution -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0 2 4 6 8 10 12 14 16 18 20  -   ( 0 ) [ s , f  (0)]= [-1.5357, 2774.813008] [ c , f  (0)]= [-3.0993, 0.511846] lower branch solution upper branch solution MOHAMMAD M. RAHMAN AND IOAN POP 56 solution f will grow enormously with time when   as can be seen from Eq. (22). Therefore, the lower branch solution is unstable and not physically realizable. Figure 3 shows that values of (0) for the upper branch solution increase with the increase of  whereas they decrease with the increase of  . As  approaches to s  , the values of (0) increase very rapidly. At s   , this value is (0) 2774.813  . In the figure we have truncated the lower branch solution as in this scale upper branch solution is not visible properly. Table 3. Smallest eigenvalues  for different values of a , b when 0    (Casson fluid) and   (Newtonian fluid) at 2.5s  , 0.5   , 0Ec  , Pr 1 . a b  Smallest eigenvalues  Upper branch, c  Lower branch, s  0 1 1 1 1 1 -1 -1 0 -1 -1 -1 2 2 2 10 100  1.5528 1.5880 1.8004 1.6276 1.6342 1.9308 -0.1262 -0.1316 -0.1443 -0.1973 -0.2169 -0.2193 In Figures 4 and 5, respectively, the effects of the suction parameter s on the reduced skin friction coefficient and Nusselt number are displayed. These figures confirm that dual solutions can be found when 0 s s s  , one solution when 0 c s s s s   and no solution when c s s . The critical 2.1001 s s  belongs to the lower branch solution while 1.6993 c s  belongs to the upper branch solution when other parameters values 1,   1,a  1,b   2  , 1Ec  , and Pr 1 are fixed. The critical 1.6993 c s  is lower than the value 2 calculated by Miklavčič and Wang [36] for a steady viscous (Newtonian) fluid flow over a permeable linearly shrinking surface. Thus, for a Casson fluid, flows over a shrinking surface a minimum fluid withdrawal from the boundary layer compared to the viscous fluid will produce dual solutions. The values of (0)f  for the upper branch solution first increase with the increase of s up to a certain value when 1cs s s  (say), then decrease with the further increase of 1 s s . On the other hand, the values of (0)f  for the lower branch solution first decrease up to a certain value with the increase of s when 2ss s s  (say), then they increase with the further increase of 2 s s . It is worth noting that 2 1 s s . The values of (0)f  for the upper branch solution are higher than those of the lower branch solution within the solution domain 4.23 s s s  . For 4.23s  an opposite trend is observed. On the other hand values of (0) (Figure 5) for the upper branch solution increase with the increase of c s s while for the lower branch solution these values first decrease to a minimum value when 3s s s s  (say) then start to increase with the further increase of the suction parameter. It is to be mentioned that to display both the solutions branch in the same scale we have truncated the lower branch solution vertically while plotting. EFFECTS OF SECOND-ORDER SLIP AND VISCOUS DISSIPATION 57 Figure 4. Variations of (0)f  for different values of s at 1   when 1a  , 1b   , 2  , 1Ec  , and Pr 1 . Figure 5. Variations of (0) for different values of s at 1   when 1,a  1,b   , 2,  1Ec  , and Pr 1 . Figure 6. Variations of (0)f  for different values of  at 1   when 1a  , 1b   , 2.5,s  1,Ec  and Pr 1 . Figure 7. Variations of (0) for different values of  at 1   when 1a  , 1b   , 2.5s  , 1,Ec  and Pr 1 . In Figures 6 and 7 we have displayed the reduced skin friction coefficient and reduced Nusselt number for various values of the Casson parameter  when 1,   1,a  1,b   2.5,s  1Ec  , and Pr 1 . For these studied parameter values two critical values of  , one for the lower branch 1.1818 s   and the other for the upper branch 0.5620, c   are identified for the existence of the solution. These figures show that dual solutions can be found when s   . The upper branch solution is always bigger than the lower branch solution. Within the region 0.85 c    the upper solution increases quite rapidly to its maximum value 0.3333 then starts to decrease with the further increase of  . As mentioned earlier a large  i.e.   corresponds to the viscous (Newtonian) fluids. Figure 6 clearly indicates that as  becomes large the physically realizable values of (0)f  approach to their asymptotic value 0.2724. Figure 7 indicates that the asymptotic value of the reduced 1.5 2 2.5 3 3.5 4 4.5 5 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 s f   ( 0 ) [s s , f  (0)] = [2.1001, 0.263762] [s c , f  (0)] = [1.6993, 0.322851] upper branch solution lower branch solution 1.5 2 2.5 3 3.5 4 4.5 5 0 1 2 3 4 5 6 7 8 9 10 s -   ( 0 ) [s c , f  (0)] = [1.6993, 0.774767] [s s , f  (0)] = [2.1001, 1428.250813] upper branch solution lower branch solution 0 1 2 3 4 5 6 7 8 9 10 0.22 0.24 0.26 0.28 0.3 0.32 0.34  f   ( 0 ) upper branch solution lower branch solution [ c , f  (0)] = [0.562, 0.298691] [ s , f  (0)] = [1.1818, 0.240036] 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10  -   ( 0 ) upper branch solution lower branch solution [ s , -  (0)] = [1.1818, 1662.724374] [ c , -  (0)] = [0.562, 1.438864] MOHAMMAD M. RAHMAN AND IOAN POP 58 Nusselt number is 2.2421 when  becomes large. This figure clearly reveals that the rate of heat transfer in a Casson fluid is lower compared to the rate of heat transfer in a viscous (Newtonian) fluid. Table 4. Critical  and corresponding values of (0)f  for different a , b , and  when 2.5s  , 1Ec  , and Pr 1 . a b  Upper branch solution Lower branch solution c  (0)f  s  (0)f  0 -1 2 2.5 5  -2.0712 -2.4046 -3.4199 -5.3429 1.027056 1.176313 1.590792 2.271937 -1.0955 -1.2444 -1.6454 -2.3477 0.440174 0.528697 0.767231 1.155035 1 -1 2 2.5 5  -3.0993 -3.5834 -5.0176 -7.6291 1.028129 1.179225 1.599951 2.293422 -1.5357 -1.7732 -2.4330 -3.5028 0.440184 0.528737 0.772045 1.155058 1 0 2 2.5 5  -2.0002 -2.2196 -2.8136 -3.7550 1.002087 1.151879 1.577269 2.281026 -2.0001 -2.2195 -1.9169 -2.5438 0.996414 1.147576 0.771982 1.155041 Table 5. Critical  and corresponding values of (0) for different a , b , and  when 2.5s  , 1Ec  , Pr 1. a b  Upper branch solution Lower branch solution c  (0) s  (0) 0 -1 2 2.5 5  -2.0712 -2.4046 -3.4199 -5.3429 0.550963 0.426726 0.090219 -0.434987 -1.0955 -1.2444 -1.6454 -2.3478 2642.053988 3003.454371 4323.338580 6507.086524 1 -1 2 2.5 5  -3.0993 -3.5834 -5.0176 -7.6291 0.511845 0.359193 -0.024985 -0.610166 -1.5357 -1.7732 -2.4330 -3.5028 2774.813008 3652.195429 4201.941189 6579.395747 1 0 2 2.5 5  -2.0002 -2.2196 -2.8136 -3.7550 0.018269 -0.193895 -0.698576 -1.568783 -2.0001 -2.2195 -1.9169 -2.5438 -0.026184 -0.242078 3387.651463 5567.562933 In Tables 4 and 5 we have investigated in details the effects of the first order slip parameter a , second order slip parameter b on the critical  , reduced skin friction coefficient and local Nusselt number for different values of the Casson parameter  when 2.5s  , 1Ec  , and Pr 1 are fixed. From Table 4 we notice that the critical | | c  for the upper branch solution as well as | | s  for the lower branch solution increases with the increase of the Casson parameter  , first order slip parameter a , and absolute value of the second order slip parameter | |b . We also notice that values of (0)f  increase with the increase of a as well as | |b for a fixed value of the Casson parameter  . Thus, the presence of first and second order slip parameters broadens the solution space before the boundary layer is going to separate. Table 5 shows that the critical | | c  increases with the increase of the Casson parameter  , first order slip parameter a , and second order slip parameter | |b . The EFFECTS OF SECOND-ORDER SLIP AND VISCOUS DISSIPATION 59 local Nusselt number ( (0) ) decreases with the increase of  and a , while it increases with the increase of | |b . For large  , the local Nusselt number becomes negative which indicates that heat transfer takes place from the fluid to the surface. Figure 8. Velocity profiles ( )f  for different values of  when 1   , 1a  , 1b   , 2.5s  , 1Ec  , and Pr 1 . Figure 9. Temperature profiles ( )  for different values of  when 1   , 1a  , 1b   , 2.5s  1,Ec  and Pr 1 . Figure 10. Velocity profiles ( )f  for different values of s when 1   , 1a  , 1b   , 2  , 1Ec  , and Pr 1 . Figure 11. Temperature profiles ( )  for different values of s when 1   , 1a  , 1b   , 2  , 1Ec  , and Pr 1 . The variations of the dual solutions in terms of dimensionless velocity '( )f  and temperature ( )  profiles of a Casson fluid for different values of the Casson parameter , suction parameter ,s second order slip parameter b , first order slip parameter ,a and Eckert number Ec are demonstrated in Figures 8 to 16 for a shrinking surface ( 1   ) when for Pr 1 . It is seen from Figure 8 that the hydrodynamic boundary layer thickness for the upper branch solution is always thinner than that of the lower solution branch. This figure also 0 1 2 3 4 5 6 7 8 9 10 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0  f  (  )  = 1.2, 2, 10, 50 upper branch solution lower branch solution  = 1.2, 2, 10, 50 0 1 2 3 4 5 6 7 8 9 10 -0.5 0 0.5 1   (  )  = 2, 10, 50  = 2, 10, 50 lower branch solution upper branch solution 0 1 2 3 4 5 6 7 8 9 10 -1 -0.8 -0.6 -0.4 -0.2 0  f  (  ) upper branch solution lower branch solution s = 2.5, 3, 3.5 s = 2.5, 3, 3.5 0 1 2 3 4 5 6 7 8 9 10 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1   (  ) upper branch solution lower branch solution s = 2.5, 3, 3.5 MOHAMMAD M. RAHMAN AND IOAN POP 60 reveals that the velocity profile ( )f  increases with the increasing values of the Casson parameter  for the upper branch solution. For smaller values of  the rate of increase in velocity is quite rapid compared to the rate of increase for large values of  . For large values of   , the fluid behaves like a viscous (Newtonian) fluid. On the other hand the temperature of the Casson fluid and hence the thickness of the thermal boundary layer decreases for the upper branch solution with the increase of  (see Figure 9). An opposite trend is observed for the lower branch solution. These figures clearly demonstrate that the thicknesses of the hydrodynamic and thermal boundary layers are thinner for a viscous (Newtonian) fluid compared to the Casson fluid. Figure 10 reveals that the velocity profile ( )f  increases, and hence the hydrodynamic boundary layer thickness decreases, with the increasing values of the suction parameter s for the upper branch solution. A reverse effect of s on the velocity profiles is observed for the lower branch solution. Figure 11 shows that the temperature profiles (𝜂) in the vicinity of the surface decrease with the increase of s either for an upper solution branch or for a lower solution branch. The thickness of the thermal boundary layer for the upper solution branch is lower than the corresponding thickness of the lower solution branch. Figure 12. Velocity profiles ( )f  for different values of b when 1   , 1a  , 2.5s  , 2  , 1Ec  , and Pr 1 . Figure 13. Temperature profiles ( )  for different values of b when 1   , 1a  , 2.5s  , 2  , 1Ec  , and Pr 1 . The dimensionless velocity '( )f  and temperature ( )  profiles of a Casson fluid due to the second order slip parameter b and first order slip parameter a are displayed in Figures 12-15, when 2  , 2.5s  , 1   , Pr 1 . It is seen that the velocity distributions ( )f  (Figures 12 and 14) of the Casson fluids within the boundary layer increase with the increase of | |b as well as with the increase of a for the upper branch solution, whereas for the lower solution branch these distributions decrease except very close to the surface of the shrinking sheet where they increase with the increase of the parameter | |b and also with the increase of a . The thicknesses of the thermal boundary layer decreases with the increase of the parameter | |b and also with ,a as can be seen from Figures 13 and 15. Thus, the second and first order slip parameters reduce the thicknesses of the hydrodynamic and thermal boundary layers. We notice that the effect of a is more pronounced than the effect of b on the velocity and temperature fields. Therefore, first and second order slips play an important role in modeling boundary layer flows with Casson fluids over a shrinking surface. Finally, the effect of the Eckert number Ec on the temperature distribution is depicted in Figure 16 when 1   , 2  , 1a  , 1b   , 2.5s  , and Pr 1 . For these studied parameter values we notice that temperature distribution within the boundary layer for the upper branch solution increases with the increase of Ec . An increasing Ec generates more frictional heating which in turn induces the flow rate to increase. An opposite effect of Ec on the temperature distribution is observed for the lower branch solution. 0 1 2 3 4 5 6 7 8 9 10 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0  f  (  ) upper branch solution lower branch solution b = 0, -0.5, -1 b =0, -0.5, -1 0 1 2 3 4 5 6 7 8 9 10 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1   (  ) upper branch solution lower branch solution b = 0, -0.5, -1 b = 0, -0.5, -1 EFFECTS OF SECOND-ORDER SLIP AND VISCOUS DISSIPATION 61 Figure 14. Velocity profiles ( )f  for different values of a when 1   , 1b   , 2.5s  , 2,  1Ec  , and Pr 1 . Figure 15. Temperature profiles ( )  for different values of a when 1   , 1b   , 2.5s  , 2,  1Ec  , and Pr 1 . Figure 16. Temperature profiles ( )  for different values of Ec when 1   , 1a  , 1b   , 2.5s  , 2  , and Pr 1 . 6. Conclusion In this paper we investigate the steady forced convective boundary layer flow and heat transfer characteristics of a Casson fluid over a permeable shrinking surface with variable temperature in the presence of second-order slip at the interface. A numerical simulation is carried out to investigate the existence of the dual solutions. The critical shrinking, suction and Casson parameters have been identified for the existence of the dual solutions. Following our numerical computations it is concluded that dual solutions exist only when 0 s    , s s s and s   for fixed values of the other parameters where the subscript s stands for the solution corresponding to the lower branch. The upper branch solution is found to be stable and hence physically 0 1 2 3 4 5 6 7 8 9 10 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0  f  (  ) a = 0, 1, 5 upper branch solution lower branch solution a = 0, 1, 5 0 1 2 3 4 5 6 7 8 9 10 -2.5 -2 -1.5 -1 -0.5 0 0.5 1   (  ) lower branch solution upper branch solution a = 0, 1, 5 a = 0, 1, 5 0 1 2 3 4 5 6 7 8 9 10 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1   (  ) Ec = 0.1, 1, 3 Ec = 0.1, 1, 3 upper branch solution lower branch solution MOHAMMAD M. RAHMAN AND IOAN POP 62 acceptable, while the lower branch is unstable and hence physically not realizable. The critical | | c  corresponding to the stable solution increases with the increase of the Casson parameter , first order slip parameter ,a and magnitude of the second order slip parameter | | .b The rate of heat transfer in a Casson fluid is lower than the rate of heat transfer in a viscous (Newtonian) fluid. The thicknesses of the hydrodynamic and thermal boundary layers are thicker for a Casson fluid than a viscous (Newtonian) fluid. The presence of first and second order slip parameters broaden the solution space before the boundary layer is given a chance to separate. An extension of this work is underway exploring the heat transfer augmentation considering Casson fluid- based nanofluid. 7. Acknowledgement We would like to express our thanks to the very competent reviewers for their valuable comments and suggestions. M.M. Rahman is grateful to The Research Council (TRC) of Oman for funding under the Open Research Grant Program: ORG/SQU/CBS/14/007 and Dr. N. Al-Salti, Department of Mathematics and Statistics, Sultan Qaboos University for translating the English Abstract into Arabic. References 1. Sparrow, E.M. and Abraham, J.P. Universal solutions for the streamwise variation of the temperature of a moving sheet in the presence of a moving fluid. Int. J. Heat Mass Trans., 2005, 48, 3047-3056. 2. Abraham, J.P. and Sparrow, E.M. Friction drag resulting from the simultaneous imposed motions of a free stream and its bounding surface. Int. J. Heat Fluid Flow, 2005, 26, 289-295. 3. Crane, L.J. Flow past a stretching plate. J. Appl. Math. Physics (ZAMP), 1970, 21, 645–647. 4. Skelland, A.H.P. Non-Newtonian Flow and Heat Transfer, John Wiley & Sons, New York, 1966. 5. Denn, M.M. Boundary-layer flows for a class of elastic fluids. Chem. Eng. Sci., 1967, 22, 395–405. 6. Rajagopal, K.R., Gupta, A.S. and Wineman, A.S. On a boundary layer theory for non-Newtonian fluids. Appl. Sci. Eng. Lett., 1980, 18, 875–883. 7. Bird, R.B., Armstrong, R.C. and Hassager, O. Dynamics of Polymer Liquids (2 nd edition), Wiley, New York, 1987. 8. Slattery, J.C. Advanced Transport Phenomena, Cambridge University Press, Cambridge, 1999. 9. Rivlin, R.S. and Ericksen, J.L. Stress deformation relations for isotropic materials. J. Rational Mech. Anal., 1955, 4, 323–425. 10. Oldroyd, J.G. On the formulation of rheological equations of state. Proc. R. Soc. London A, 1950, 200, 523–541. 11. Hayat, T., Shafiq, A. and Alsaedi A. Effect of Joule heating and thermal radiation in flow of third grade fluid over radiative surface. PLoS ONE, 2015, 9(1), e83153, doi:10.1371/journal.pone.0083153. 12. Fredrickson, A.G. Principles and Applications of Rheology. Prentice-Hall, Englewood Cliffs, New Jersey, 1964. 13. Eldabe, N.T.M. and Salwa, M.G.E. Pulsatile magnetohydrodynamic viscoelastic flow through a channel bounded by two permeable parallel plates. J. Phys. Soc. Japan,1995, 64, 4163-4174. 14. Boyd, J., Buick, J.M. and Green, S. Analysis of the Casson and Carreau-Yasuda non-Newtonian blood models in steady and oscillatory flow using the lattice Boltzmann method. Phys. Fluids, 2007, 19, 93-103. 15. Mernone, A.V., Mazumdar, J.N. and Lucas, S.K. A mathematical study of peristaltic transport of a Casson fluid. Math. Comp. Model., 2002, 35, 895–912. 16. Mustafa, M., Hayat, T., Pop, I. and Aziz, A. Unsteady boundary layer flow of a Casson fluid due to an impulsively started moving flat plate. Heat Transfer Asian Res., 2011, 40, 563–576. 17. Pramanik, S. Casson fluid flow and heat transfer past an exponentially porous stretching surfcae in presence of thermal radiation. Ain Shams Eng. J., 2014, 5, 205-212. 18. Nakamura, M. and Sawada, T. Numerical study on the flow of a non-Newtonian fluid through an axisymmetric stenosis. ASME J. Biomechanical Eng., 1988, 110, 137–143. 19. Wu, L. A slip model for rarefied gas flows at arbitrary Knudsen number. Appl. Phys. Lett., 2008, 93, 253103 (3 pages). 20. Fang, T., Yao, S., Zhang, J. and Aziz, A. Viscous flow over a shrinking sheet with a second order slip flow model. Commun. Nonlinear Sci. Numer. Simulat., 2010, 15, 1831–1842. 21. Mahmood, T., Shah, S.M. and Abbas, G. Magnetohydrodynamic viscous flow over a shrinking sheet with second order slip flow model. Heat Trsnf. Res., 2015, doi:10.1615/HeatTransRes.2015007512. 22. Fukui, S. and Kaneko, R. A database for interpolation of Poiseuille flow rates for high Knudsen number lubrication problems. ASME J. Tribol., 1990, 112, 78–83. 23. Weidman, P.D., Kubittschek, D.G. and Davis, A.M.J. The effect of transpiration on self-similar boundary layer flow over moving surfaces. Int. J. Eng. Sci., 2006, 44, 730-737. EFFECTS OF SECOND-ORDER SLIP AND VISCOUS DISSIPATION 63 24. Roşca, A.V. and Pop, I. Flow and heat transfer over a vertical permeable stretching/shrinking sheet with a second order slip. Int. J. Heat Mass Trans., 2013, 60, 355–364. 25. Roşca, N.C. and Pop, I. Mixed convection stagnation point flow past a vertical flat plate with a second order slip: heat flux case. Int. J. Heat Mass Trans., 2013, 65, 102-109. 26. Rahman, M.M., Roşca, A.V. and Pop, I. Boundary layer flow of a nanofluid past a permeable exponentially shrinking/stretching surface with second order slip using Buongiorno’s model. Int. J. Heat Mass Trans., 2014, 77, 1133-1143. 27. Rahman, M.M. and Pop, I. Mixed convection boundary layer stagnation-point flow of a Jeffrey fluid past a permeable vertical flat plate. Zeitschrift fuer Naturforschung A (ZNA), 2014, 69(12), 687-696. 28. Rahman, M.M., Roşca, A.V. and Pop, I. Boundary layer flow of a nanofluid past a permeable exponentially shrinking surface with convective boundary condition using Buongiorno’s model. Int. J. Num. Methods Heat & Fluid Flow, 2015, 25(2), 299-319. 29. Harris, S.D., Ingham, D.B. and Pop, I. Mixed convection boundary-layer flow near the stagnation point on a vertical surface in a porous medium: Brinkman model with slip. Trans. Porous Media, 2009, 77, 267-285. 30. Shampine, L.F., Gladwell, I. and Thompson, S. Solving ODEs with MATLAB. Cambridge University Press, Cambridge, 2003. 31. Shampine, L.F., Reichelt, M.W. and Kierzenka, J. Solving boundary value problems for ordinary differential equations in MATLAB with bvp4c, 2010. . 32. Wang, C.Y. Flow due to a stretching boundary with partial slip-an exact solution of the Navier-Stokes equations. Chem. Eng. Sci., 2002, 57, 3745-3747. 33. Wang, C.Y. Analysis of viscous flow due to a stretching sheet with surface slip and suction. Nonlinear Anal. Real World Appl., 2009, 10(1), 375-380. 34. Anderson, H.I. Slip flow past a stretching surface. Acta Mech., 2002, 158, 121-125. 35. Sahoo, B. and Do, Y. Effects of slip on sheet-driven flow and heat transfer of third grade fluid past a stretching sheet. Int. Commun. Heat Mass Transfer, 2010, 37, 1064-1071. 36. Miklavčič, M. and Wang, C.Y. Viscous flow due to a shrinking sheet. Q. Appl. Math., 2006, 64, 283-290. Received 18 March 2015 Accepted 10 November 2015