Effects of Second-Order Slip and Viscous Dissipation on the Analysis of the Boundary Layer Flow and Heat Transfer Characteristics of a Casson Fluid

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.


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 twodimensional 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.

Basic Equations
Consider the steady two-dimensional flow of a Casson fluid past a permeable shrinking sheet with the velocity 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 Txwhile the constant temperature of the ambient (inviscid) fluid is , T  where ) ( (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, or, (see Nakamura and Sawada [18]) 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, ) 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 and  is the momentum accommodation coefficient with 1 0    .Based on the definition of l , it is seen that for any given value of n K , we have . 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: where  is the stream function, which is defined in the classical form as where prime denotes differentiation with respect to  .Substituting (7) into Eqs.( 4) and ( 5), we obtain the following ordinary (similarity) equations subject to the boundary conditions where a is the first order velocity slip parameter with , 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 x Nu , which can be easily shown to be given by where Re ( ) /

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 subject to the initial and boundary conditions 0 0 slip 0 : , 0, for any , 0 : , ( ) ( ), ( ) at 0 0, as where t denotes time.

Numerical Technique
Following Rahman et al. [26][27][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 ( , ) subject to the two-point boundary conditions ( ( ), ( )) 0 bc x a x b  .(30) The approximate solution () t  is a continuous function that is a cubic polynomial on each subinterval  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 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  , 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) 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    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.

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.     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. 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 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    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

s 
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 ss  while for the lower branch solution these values first decrease to a minimum value when 3 s 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.
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,   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  In Tables 4 and 5      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.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.

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 ss  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

IFigure 1 .
Figure 1.The geometry of the problem of a shrinking surface.

,
where y is the coordinate measured in the normal direction to the surface and the coordinate 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 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, problem (19)-(21), we write (see Weidman et al.[23] or Rahman et al.[26][27][28]) Introducing(22) into Equations.(

.
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 For the present problem, we relax the condition that 0 ( ) 0

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  40-60 for the lower branch (second) solution.
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 0 s  and second order slip parameter 0 b  , equation (10

for 1 a
 and several values of s and b when 1

Figure 2
Figure 2. Variations of (0) f  for different values

Figure 3
Figure 3. Variations of (0)    for different values

.
In the figure we have truncated the lower branch solution as in scale upper branch solution is not visible properly.
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 1 c s s s (say), then decrease with the further increase of 1 ss  .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 2 s s s s  (say), then they increase with the further increase of 2 ss  .It is worth noting that 21 ss  .The values of (0) f  for the upper branch solution are higher than those of the lower branch solution within the solution domain 4

Figure 4 .
Figure 4. Variations of (0) f  for different values of

Figure 5
Figure 5. Variations of (0)    for different values

and Pr 1 
. For these studied parameter values two critical values of  , one for the lower branch 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.
 (0)] = [0.562,1.438864]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.

1 
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 2are 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 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.

Table 1 .
Comparison of results

Table 2 .
Comparison of Table 3 presents the smallest eigenvalues  for the upper branch solution at several values of , a b and .

Table 3 .
Smallest eigenvalues  for different values of

Table 4 .
Critical  and corresponding values of (0) f 

Table 5 .
Critical  and corresponding values of (0) 