Hydrodynamic flow of non-Newtonian power-law fluid past a moving wedge or a stretching sheet: a unified computational approach

A unified mathematical equation that combines two different boundary-layer flows of viscous and incompressible Ostwald-de Waele fluid is derived and studied. The motion of the mainstream and the wedge is approximated in the power-law manner, i.e, in terms of the power of the distance from the leading boundary-layer edge. It is also considered that the wedge can move in the same and opposite direction to that of the mainstream. The governing partial differential equations are transformed into the nonlinear ordinary differential equation using a new set of similarity variables. This transformed equation subjected to the boundary conditions describing the flow is then solved using the Chebyshev collocation method. Further, these numerical results are then validated by determining the flow behaviour at far-field by performing asymptotics. The velocity ratio parameter effectively captures and distinguishes two boundary-layer flows. The boundary layer thickness for shear-thinning fluid is thinner compared to corresponding parameters for shear-thickening fluids and is markedly separated by the Newtonian fluid. Further, the boundary-layer flow of the non-Newtonian fluid predicts an infinite viscosity for shear-thinning fluid quite away from the surface. The hydrodynamics of the obtained results is explained thoroughly.

www.nature.com/scientificreports www.nature.com/scientificreports/ equations look similar but their description of motion is entirely different. Denoting U x ( ) and U x ( ) w as the velocity of the mainstream and the wedge surface respectively, let ε be the ratio of the wedge velocity to the mainstream velocity. The parameter ε mathematically simplifies and unifies the above mentioned two boundary-layer flow problems for ε = 0 and ε = 1 (see section 2 for details). In other words, the mathematical equation for ε = 0 that represents the boundary layer flow due to a stretching surface whose solution is substantially different from the boundary-layer flow over a constant wedge (ε = 1). In this line of studies, Afzal 31 derived a new type of boundary-layer flow problem by modifying the similarity transformations that produce ε into it and solved numerically. Kudenatti 32 revisited the same problem and solved it by uniformly valid convergent series solution in the range of ε ∈ . [0, 0 5]. Recently, Karkera et al. 33 extended it by including the magnetic field and attempted the Haar wavelet-based numerical technique and compared their results with existing literature for various ε values. In the present problem, a unified mathematical model is again derived for the hydrodynamic boundary-layer flow of a power-law fluid (Ostwald-De Waele) over a wedge, since the flow of non-Newtonian fluid has ample applications in chemical industries.
There are numerous applications of power-law fluids in various manufacturing industries. Sewage sludge, china clay, oil-water emulsion, cosmetics, paints, synthetic lubricants, biological fluids, jam, jellies, etc are some of the non-Newtonian fluids (Schowalter 34 ; Bird et al. 35 , Andresson & Irgens 36 ; Shenoy & Mashelkar 37 , etc). For these fluids, the viscosity or apparent viscosity essentially depends on the shearing rate. Acrivos et al. 38 theoretically analyzed the power-law fluid flow and heat transfer over a surface considering various external surfaces. Mitschka & Ulbrecht 39 and Andersson et al. 40 have numerically presented the basic solutions for the shear-thickening and shear-thinning power-law fluids. However, the solutions shown by these authors did not match the mainstream properly. This issue is later fixed by Denier & Hewitt 41 by modifying the similarity solutions of the boundary-layer flow for shear-thinning and shear-thickening fluids. Ishak et al. 42 considered the boundary-layer flow of power-law fluid over a moving wedge for which the model shows non-unique solution structures which are similar to the Newtonian case. The linear stability of Griffith et al. 43 shows that the stationary spiral instabilities observed experimentally can be exactly described by shear-thinning fluids. Griffiths 44 has discussed the stability of shear-thinning boundary-layer flow of the power-law fluids in which the base flow is obtained using the similarity transformations of the Prandtl boundary layer equations for non-Newtonian fluid. He further discusses that the viscosity of the power-law fluid remains unbounded for shear-thinning fluid which is unphysical. Denier & Dabrowski 45 have identified a non-unique solution to the shear-thickening and shear-thinning fluid flow in the boundary-layer. They further have shown asymptotically that the decaying of shear-thinning solutions into the mainstream is strongly algebraic.
In the present study, an attempt is made to predict the steady boundary-layer flow behaviour of a non-Newtonian power-law fluid over a stretching surface or a moving wedge whose speed is approximated in terms of the power of the distance from the leading edge of the boundary-layer. This work may be regarded as an extension of the flow problem encountered by Kudenatti 32 to a class of non-Newtonian fluid known as Ostwald-de-Waele fluid. Since the resultant equation is highly nonlinear, we employ the Chebyshev collocation method(CCM) for the computation of the pertinent results. The CCM is a powerful technique to predict flow behaviour in the boundary-layer. There are very few papers that use CCM in the context of boundary-layer, however, for nonlinear or mild nonlinear see 46 .
Organization of the paper is as follows. We discuss the formulation of the problem in section 2 and the self-similar solutions are computed with suitable similarity transformations. The methodology applied to solve the flow problem is discussed in detail in section 3. To assess the nature of these obtained numerical solutions, we predict the flow behaviour of the system asymptotically from a large distance away from the origin of the boundary layer in section 4 known as far-field behaviour. Section 5 enunciate the results and the corresponding hydrodynamics is discussed in detail and a brief conclusion of the study is made in section 6.

Formulation
We study the laminar flow of a steady, incompressible non-Newtonian power-law fluid over a semi-infinite moving wedge or due to a stretching sheet in a two-dimensional space. Let the streamwise coordinate along the wall be x and normal to the wall be z (refer Fig. 1 for detailed geometric representation). The flow is governed by the two-dimensional Navier-Stokes equation accompanied by continuity equation where v = u v ( , ) is the velocity vector, ρ is the density of the fluid, p is the pressure, τ is the deviatoric stress tensor. The stress tensor for the power-law fluid is depicted by the relation where γ  is the second invariant of the strain-rate tensor and the constitutive viscosity relation is described by the Ostwald-de Waele power-law model (Bird et al. 35 ) Here K and m respectively denote the consistency index and the degree of non-Newtonian behavior of the fluid (Denier & Dabrowski 45 ). The equations governing the model are non-dimensionalized with appropriate dimensionless variables where L, U are appropriate scale factors for reference quantities to length and mainstream velocity respectively and δ denotes boundary-layer thickness. The non-Newtonian viscosity µ is referred to as an apparent viscosity µ app , wherein Further, the dimensionless apparent Reynolds number for the non-Newtonian fluid exhibits the form m m 2 which reduces to the Newtonian fluid when = m 1. To this end, the two-dimensional boundary layer flow of non-Newtonian fluid is studied in which a single mathematical equation unifies two types of boundary-layer flows: over a moving wedge of an included angle πβ and due to nonlinear stretching of the surface. In the former case, for the flow of mainstream fluid = U U x ( ) at high Reynolds number the viscosity effects are confined to a small region near the immediate vicinity of the wall, known as the boundary-layer region (Yuan 47 ;Kudenatti 5 ). In the latter situation, the boundary-layer forms in a still fluid ( = U x ( ) 0) due to stretching of the surface with large Reynolds number (Acheson 48 ). These two flow problems have been mathematically combined into a single equation that is derived for the first time in the literature. Substituting (8) in (6) and making use of the fact that Re is large and using the other boundary layer approximations and retrieving back to original variables, we get www.nature.com/scientificreports www.nature.com/scientificreports/ ρ ρ The conditions appropriately describing the flow are The boundary conditions (11) infer that the velocity of the fluid over the wedge surface asymptotically approaches the freestream far away from the wedge. Here U x ( ) is the mainstream velocity and U x ( ) w is the velocity with which the wedge or stretching sheet moving in the boundary layer and approximating these velocities in a power-law manner n w n 0 where ∞ U U , 0 and n are arbitrary constants. The above forms exhibit a self-similar solutions in the boundary-layer in the Ostwald-de Waele fluid (Schowalter 34 ). Further, from Bernoulli's theorem for an incompressible flow of non-Newtonian fluid over a moving wedge, we have with = u U on the edge of the boundary layer for moving mainstream (13) Consequently, the momentum Eq. (10b) becomes, in general (14) and otherwise (14) holds good for moving fluid. Accordingly, we define a single unified set of similarity transformations r q q where A, q are unknowns to be determined and U r is a suitable reference velocity defined as is the normalized stream function and η is the similarity variable. Defining the stream function such that it satisfies continuity Eq. (10a) identically to obtain a highly nonlinear third-order differential equation is the generalized pressure gradient parameter, ε = = is the velocity ratio parameter and the boundary conditions enclosing the system are We can enunciate the following cases depicted in Table 1 from (18) and (19) for different parameters affiliated with the system. From the table, it is clear that the governing system (18) and (19) unifies many physical models which are presented in the literature.
In this paper, an attempt has been made to study the non-Newtonian fluid effects on the two-dimensional boundary layer flow for different β and ε. We note that ε > 1 and ε < 0 cases correspond respectively to wedge is moving faster and slower than mainstream velocity in opposite direction, while ε ∈ .
. (0 0, 0 5) cases are analogous as above but in the same direction. At this stage, it is also to note that for > n 0 and < n 0 cases correspond to accelerated and decelerated pressure gradient in the boundary-layer while = n 0 is the Blasius flow of non-Newtonian fluid. The flow of Newtonian fluid ( = m 1) clearly demarcates the shear-thinning fluids ( < m 1) from the shear-thickening fluids ( > m 1). We employ the Chebyshev collocation method(CCM) for solution of the system (18) and (19) in which CCM effectively captures the fractional powers ( − m 1), the details of CCM are given below.

Computational procedure
In recent years, the Chebyshev collocation methods are gaining wide popularity since the methods are robust and provide solutions which complements precisely the physical nature of the flow. While the CCM is well known for linear or mildly nonlinear fluid mechanics problems, the application of CCM to the problem in question is non-standard. We stress that the problems involving non-Newtonian fluid, we have not observed how to apply CCM. Thus, briefing the details of CCM is definitely of interest in its own right. We discuss the methodology below.
A standard i th order nonlinear ordinary differential equation can be written in the form where P, Q,  and V are known continuous functions of η in a b [ , ] and is supplied with necessary boundary conditions. With the help of a suitable linear transformation the interval a b [ , ] can be transformed to − [ 1,1] which is the domain of Chebyshev polynomials η T ( ) n (Boyd 49 ). Then η y( ) and its derivatives have a truncated Chebyshev series of the form where a n and a n i ( ) are unknown Chebyshev coefficients of η y( ) and its derivatives. Let η j denote the Chebyshev collocation points defined as where N is a positive integer chosen at which the free stream condition is satisfied. The above Eqs. (21) and (22) have a matrix representation with www.nature.com/scientificreports www.nature.com/scientificreports/ where ′ denotes a transpose of the quantities. Further, it readily follows that Using the relation by coefficient matrix of η y( ) and its derivative ( (27) can be written as and these are evaluated at Chebyshev collocation points. Further simplifying using (29), . Consequently, the matrix form of (20) can be written as The above equation represents a system of + N ( 1) nonlinear algebraic equations involving unknown Chebyshev coefficients A given by (26). Further, incorporating the boundary conditions, we obtain the modified augmented matrix which on inversion gives (26), the required solution can be computed from (21). On the other hand, to solve the flow equation of the model considered here (18) and (19), we first transform the semi-infinite flow domain ∞ [0, ) to − [ 1,1] by using the following transformations where η ∞ is a large infinity value at which the boundary condition is met. Substitution of (35) in (18) and (19) gives nonlinear ordinary differential equation in the interval − [ 1,1] for which the solution can be written in terms of the truncated Chebyshev series as www.nature.com/scientificreports www.nature.com/scientificreports/ where X is the unknown coefficient vector of the Chebyshev polynomials. Following (22), (29), (30) and (31), we have the final matrix form as 1 when ≠ m 1. We can overcome the latter concern by performing simple computational linear algebra. A close examination reveals that the matrix which is obtained as a product form of TM X 2 2 2 is diagonalizable. Hence, it can be written in terms of its eigenvalues D and eigenvec- in all the simulations and our scientific code essentially performs the diagonalization. Hence we can write the augmented matrix as = WX S (38) where


We now incorporate the boundary conditions (19) by using (29).
We modify Eq. (38) by plugging the row matrices (39) at first, second and last row of (38) respectively. Therefore, we obtain the modified augmented matrix of the form Since X is unknown, the left hand side of (37) represents the augmented form of non-linear algebraic equations and this can be handled by taking an initial approximation to X (in the present computation X is taken as a zero vector of order + N ( 1, 1) chosen at the first phase of computation) and updating it at each stage of calculation for sufficient N until a desired tolerance of − 10 6 is attained. We observed that, the number of terms required to achieve the desired results were different for different m. Suppose = m 1 we required a maximum of 30 terms and for < m 1 say = .
. m 0 6, 0 8 an approximate of 27 terms were sufficient to get the desired accuracy whereas for > m 1 say = .
. m 1 2, 1 4 number of terms required were 40 (refer to Fig. 2 and Table 2 for details). All the results provided below are tested for their convergence criteria. Solving Eq. (40) for X, and hence retrieving to original domain by using (35) in (36) we can further compute the derivatives of (36) to obtain the velocity profiles and skin-friction coefficient associated with the Ostwald-de Waele flow. Figure 2 also emphasizes the convergence of the CCM as the computations are progressed with increasing the number of terms. At Fig. 2a, the curve is not sufficiently smooth, although conditions are clearly satisfied, but Fig. 2d simulates the required smooth profile. Table 2 also gives the corresponding wall shear stress values at each N values. From the table, it is observed that for = .
m 0 6, the convergence occurs much before 40 terms in the CCM, but to make table consistent it is given upto 40. On the other hand, the CCM is also tested by fixing the number of terms equal to 50 and decreasing the tolerance value to − 10 8 . The results are indistinguishable in both cases, hence, we continued all our simulations with the former analysis.

Large η asymptotics
In order to validate the presented numerical results achieved using CCM, we also perform asymptotic analysis for the system (18) and (19) when η  1 with the objective of having an exact solution in terms of confluent hypergeometric functions. We follow the similar procedure given by (Sachdev et where η 0 is an arbitrary constant, η S( ) and its derivatives are assumed small. Note that η = 0 0 from (19). The above form (41) is valid for all m, n and ε and η S( ) tends to zero away from the wedge surface. Substitution of (41) in (18) and (19) and upon linearizing, we have with the conditions Using the transformations   , ⋅ ⋅ ⋅ U( , , ) and ⋅ ⋅ ⋅ L( , , ) are the confluent hypergeometric function of second kind and Laguerre function respectively, and c 1 and c 2 are arbitrary constants (Abramowitz & Stegun 51 and Andrews 52 ). The special functions U and L can be further transformed into the confluent hypergeometric function of first kind using (  (47), (46) can be rewritten in original variables (44) as  (48) gives To eliminate c 2 , we utilize another relationˆˆˆˆˆ= In view of (51) for η = →∞ ( 1) for zero-pressure gradient = n 0 and ε = 1. This form certainly predicts the far-field behavior, but with larger domain. Nevertheless, the solution given by (53) is valid for all m, n and ε and effectively captures (54) as a special case. The solution for η S( ) and hence the velocity profile is given by (53). Some of the velocity profiles φ η ′( ) obtained from (53) are shown in Fig. 3 for different ε and m values. From Fig. 3, it is noticed that, all the velocity shapes are decaying to the mainstream flow asymptotically. This validates the assumption (41) that η S( ) in (52) tends to zero as η → ∞. The solution (53) takes smaller domain to satisfy the end condition which assures the exponential-type decay of the solutions. The expression (54) given by (Denier & Dabrowski 45 ) is valid for < m 1 and decays to the mainstream algebraically. However, the solution (53) is valid even for > m 1, = m 1, < m 1 as is clearly seen in Fig. 3 and for any ε. For ε = 1, it is further observed that the boundary layer thickness is found to be decreasing for increasing m value. This is in agreement with those given by Agarwal et al. 53 who have studied the boundary layer flow of non-Newtonian fluid over a thin needle.

Numerical results and discussion
The velocity ratio parameter ε = + ∞ ∞ ( ) mathematically combines two different boundary layer problems. Two extreme cases ε = 0 and ε = 1 correspond to boundary layer flow due to stretching surface and over a constant wedge. For other values of ε, boundary-layer flow over a moving wedge admits a class of self-similar solutions even for non-Newtonian fluid. When ε > 1 and ε < 0, the mainstream flow has a slower and faster velocity than the wedge and both are moving in opposite directions. However, when ε < 1 and ε > 0, both wedge and mainstream move in the same direction. While for ε = . Further, these numerical simulations obtained by CCM can be matched comparatively with results obtained by Kudenatti 32 in case of Newtonian fluid ( = m 1 in (18)). Table 3 compares various results obtained by the CCM www.nature.com/scientificreports www.nature.com/scientificreports/ with the analytical solution given by Kudenatti 32 in terms of the wall shear stress φ″(0) for various ε, β and for = m 1 (Newtonian case). An iterative scheme in CCM is used in which the velocity shapes and the wall shear stress are determined i.e, the number of terms required to produce these results is then adjusted in response to the error-tolerance fixed for the simulation. The results produced by the CCM agreed well with those given via analytical solutions in the range of ε ∈ . (0, 0 5). It is further noticed that the accuracy is found to be better for increasing ε and β. The φ″(0) value becomes zero and changes its sign at ε = .
0 5 for all β and m. We also note that the convergent solution φ″(0) obtained for certain parameters is used as an initial condition for consistency in the simulations. Table 3 also ensures that there exists a velocity profile for each value which is benign and satisfies the boundary conditions asymptotically.
We now study the broader structure of the boundary layer flow over a moving wedge for various velocity ratio parameter ε and non-Newtonian index parameter m and are given in terms of the wall shear stress φ″(0). These results are just analogous to those presented in Table 3  In other words, for increasing ε from 0, φ″(0) also increases. When ε < 1 2 , the wall shear stress values are higher for smaller n while for ε ≥ 1 2 , the trend is reversed which is observed clearly in the figure. The similar trend is observed for all m and n. Karkera et al. 33 have shown, for Newtonian fluid, similar solution structure including the double solutions when the pressure gradient parameter n is negative.
In order to realize the velocity profiles in the two-dimensional boundary layer in the power-law fluids, we have plotted the velocity shapes φ η ′( ) in Fig. 7 as a function of η for ε = . 0 2 and ε = . 0 8 and for shear-thinning and shear-thickening fluids. It is clearly seen that the velocity profiles satisfy the boundary conditions asymptotically as η → ∞. These shapes further confirm that the shear-thinning fluid < m 1 makes the boundary layer thickness thinner compared to the shear-thickening fluid > m 1. These two cases are clearly separated by the Newtonian fluid = m 1. The results are further confirmed by the asymptotic solution (53) that are given in Fig. 3 although for different value ε values. Figure 8 depicts the velocity profiles obtained by the CCM for different ε ∈ [0, 1] including two extreme cases ε = 0 and ε = 1 for shear-thinning and shear-thickening fluids. These shapes have been simulated using the CCM. The quality rheological interpretation of the solution for different ε tested and for shear thinning fluid deviates slightly from the shear-thickening fluid, which is clearly seen in Fig. 8. Around ε = .
0 5, the curves criss-cross at η ∼ . 0 5 and the criss-cross is shifting towards the wedge surface for shear-thickening fluids. Yurusoy 54 has shown for ε = 0 that exactly similar velocity shapes for shear-thinning and shear-thickening fluids wherein the velocity profiles approach the mainstream flow exponentially for ≤ m 1 and for > m 1, the shapes take larger distance to approach the mainstream. However, the CCM effectively captures these characteristics of the profiles in a smaller domain as seen in the figure. For two values of ε, some of the fluid velocity curves have been reproduced in Fig. 7 when m is varied to study the Ostwald-de Waele fluid effects on the boundary layer flow. These results are analogous to those produced in Fig. 3 which are for = .
m 0 8 and = . m 1 2 crossing the Newtonian case = m 1 at some η .  0 7 and the similar structure appears for ε = .
0 8 which is seen in the figure. However, the asymptotic solutions are given by (53) in Fig. 3 have no cross-over, this may be attributed to the exact solutions given in terms of the confluent hypergeometric functions.  www.nature.com/scientificreports www.nature.com/scientificreports/ . For a Newtonian fluid, i.e, for = m 1, µ = 1, a constant viscosity is predicted which is obvious. Two representative values of ε, say ε = 0 and ε = 1 are chosen and for a favorable pressure gradient parameter = n 1. From Fig. 9, it can be readily observed that, for the shear-thinning fluids ( = .
m 0 6 and = . m 0 8), there is quite large viscosity within the confinement of the boundary layer region. This is considered to be unphysical for the power-law fluids. An   www.nature.com/scientificreports www.nature.com/scientificreports/ enhancement in viscosity is still larger for ε = 0 compared to ε = 1 case. On the other hand, for shear-thickening fluids ( = . m 1 2 and = . m 1 4), a finite viscosity is predicted which approaches to zero as η → ∞. This clearly shows the viscosity effects are completely neglected when an inviscid flow is approached. We noticed, although not shown here, that for the increasing pressure gradient, the viscosity on the surface (µ(0)) is found to be decreasing for shear-thinning fluids. It is also noticed that the viscosity increases unboundedly away from the wedge surface as the pressure gradient β increases. While the nature of µ(0) and the viscosity has a reverse trend for the shear-thickening fluids.

Conclusions
The present work deals with the hydrodynamic flow a non-Newtonian power-law fluid (Ostwald-de Waele fluid) model past a moving wedge or a stretching sheet in a unified manner. The wedge/sheet is assumed to move/ stretched with a velocity varying as the power of the distance from the leading edge of the boundary-layer along with the freestream. With the help of suitable similarity transformations, the governing partial differential equation system was then reduced to the non-linear ordinary differential equation and then solved numerically by Chebyshev collocation technique. The application of the collocation method for the problem under consideration is first of its kind hence we propose an asymptotic solution which is in validation with these numerical results. The model was discussed for the shear-thinning and shear-thickening power-law fluids, and a comparison was made with the Newtonian solutions. Results for the skin friction coefficient, velocity profiles as well as viscosity profiles are presented for different values of the governing parameters. It is noticed that the wall-shear stress increases with an increase in ε, the velocity ratio parameter and decrease in the power-law index m. Further, the boundary-layer thickness for shear-thinning fluids is small in comparison with that of shear-thickening power-law fluids for a given set of parameters. Physically, the flow for < m 1 is always adhered to the wedge surface thereby making the flow benign whereas for > m 1, the flow is convected towards the mainstream. It is also shown that with the increase in shearing stresses, the viscosity of the shear-thinning power-law fluids < m 1