University of Birmingham Flow of a generalised Newtonian fluid due to a rotating disk

The boundary-layer ﬂow due to a rotating disk is considered for a number of generalised Newtonian ﬂuid models. In the limit of large Reynolds number the ﬂow inside the three-dimensional boundary-layer is determined via a similarity solution. Results for power-law and Bingham plastic ﬂuids agree with previous investigations. We present solutions for ﬂuids that adhere to the Carreau viscosity model. It is well known that unlike the power-law and Bingham models the Carreau model is applicable for vanishingly small, and inﬁnitely large shear rates, as such we suggest these results provide a more accurate description of non-Newtonian rotating disk ﬂow. (cid:2) 2015 The Author. Published by Elsevier B.V. ThisisanopenaccessarticleundertheCCBYlicense(http:// creativecommons.org/licenses/by/4.0/).


Introduction
The steady incompressible flow induced by the rotation of an infinite plane with uniform angular velocity is an exact solution of the Navier-Stokes equations, as was first described by von Kármán [1]. The flow is characterised by the lack of a radial pressure gradient near to the disk to balance the centrifugal forces so the fluid spirals outwards. The disk acts as a centrifugal fan, the fluid emanating from the disk being replaced by an axial flow directed back towards the surface of the disk.
Batchelor [2] showed that this type of flow is in fact just a limiting case of a whole number of flows with similarity solutions in which both the infinite plane and the fluid at infinity rotate with differing angular velocities. The corresponding limiting case, when the infinite plane is stationary and the fluid at infinity rotates at a constant angular velocity, was first described by Bödewadt [3].
A vast wealth of material exists concerning the solutions of the Newtonian rotating disk equations; the interested reader is referred to Zandbergen and Dijkstra [4]. The authors provide a thorough review of the major contributions made postdating von Kármán's seminal work.
Considerably less attention has been given to the corresponding non-Newtonian rotating disk problem. Mitschka [5] modified the von Kármán similarity solution to incorporate a power-law governing viscosity relationship. In this case the base flow is not an exact solution of the generalised Navier-Stokes equations and a boundary-layer approximation is required. Both Mitschka and Ulbrecht [6] and Andersson et al. [7] present basic flow solutions for shear-thickening and shear-thinning power-law fluids. However, the authors overlooked the importance of matching these boundary-layer solutions to an external flow. Denier and Hewitt [8] addressed this problem and presented corrected solutions for both cases, noting that the structure of the solutions is intrinsically different for shear-thickening and shear-thinning fluids.
More recently, Ahmadpour and Sadeghy [9] (subsequently referred to herein as AS) formally addressed the problem of the flow due to a rotating disk when one considers Bingham plastic fluids. Claiming to have found an exact solution to the problem, the authors are only able to present numerical solutions for specific values of the Reynolds number (Re) and dimensionless radius of the disk (r). Having not considered the boundary-layer formulation of the problem, the authors find that terms dependent on Re and r appear in the formulation of the governing base flow ODEs, and thus have the need to provide specific values for these constants during their numerical solution process.
In this study we determine steady mean flow solutions for power-law, Bingham and Carreau fluid models. The power-law results are essentially a review of the work of Denier and Hewitt [8] but are included here as they prove useful to compare with the results owing from the more complex Carreau model. By introducing the modified Bingham number used by Matsumoto et al. [10] in their film thickness investigation, we are able to determine a governing set of ODEs dependent solely on this parameter, these results are then compared to those of AS. Additionally, we present solutions for shear-thickening and shear-thinning Carreau fluids where now the flow is controlled by not one, but three dimensionless parameters. In Section 2 we formulate the problem in the general case, results are presented in Section 3 and are discussed in Section 4. Conclusions are drawn in Section 5.

Formulation
Consider the flow of a steady incompressible generalised Newtonian fluid due to a rotating disk located at z Ã ¼ 0. The disk rotates about the z Ã -axis with angular velocity X Ã . Working in a reference frame that rotates with the disk, the continuity and Cauchy momentum equations are expressed as Here u Ã ¼ ðu Ã ; v Ã ; w Ã Þ are the velocity components in cylindrical polar coordinates ðr Ã ; h; z Ã Þ, the angular velocity vector is X Ã ¼ ð0; 0; X Ã Þ, the position vector is r Ã ¼ ðr Ã ; 0; z Ã Þ, the fluid density is q Ã and p Ã is the fluid pressure. The stress tensor s Ã for incompressible generalised Newtonian fluids is given by Here m Ã is the consistency coefficient and n is the fluid index, for n > 1 the fluid is said to be shear-thickening, whilst for n < 1 the fluid is said to be shear-thinning. The Newtonian viscosity relationship is recovered when n ¼ 1; s Ã y ¼ 0 and l Ã 0 ¼ l Ã shear-rate viscosity is l Ã 1 , the zero-shear-rate viscosity is l Ã 0 and k Ã is the characteristic time constant, often referred to as the 'relaxation time'.
Assuming the flow to be axisymmetric the components of the stress tensor are where the magnitude of the rate-of-strain tensor takes the form In the rotating frame of reference this system is closed subject to the boundary conditions We non-dimensionalise the system by writing We note here that the axial coordinate and velocity component have been scaled by the boundary-layer thickness, d Ã , this is in anticipation of a boundary-layer structure arising on the rotating disk. One finds that where throughout the forthcoming analysis q ¼ n for power-law fluids and q ¼ 1 for Bingham plastic and Carreau fluids, whilst u @u @r þw @u @z À ðv þrÞ 2 r ¼À @p @r þ @ @z l @u @z þ 1 Re 2=ðqþ1Þ 2 r @ @r lr @u @r þ @ @z l @w @r À 2lu u @v @r þ w @v @z þ uv r þ 2u ¼ @ @z l @v @z þ 1 Re 2=ðqþ1Þ 1 r 2 @ @r u @w @r þ w @w @z ¼ ÀRe 2=ðqþ1Þ @p @z þ 1 r @ @r lr @u @z þ 2 @ @z l @w @z þ 1 Re 2=ðqþ1Þ 1 r @ @r where the dimensionless viscosity functions (in the yielded region when considering Bingham plastic fluids) l are defined as Power-law model À l ¼ ðlÞ nÀ1 ; ð8eÞ where B n ¼ s Ã is the Bingham number defined by Matsumoto et al. [10]. The authors experimental investigations have shown this quantity to be Oð1Þ for flows with Re ) 1. The Carreau viscosity ratio is is the dimensionless equivalent of the constant k Ã . Herê the higher order terms, Ll, that contribute to the generalised viscosity, are given in the Appendix A for completeness. We note that the expressions for B n and k can be simplified when one considers the modified Reynolds number, R ¼ rRe 1=2 , based on the boundary-layer thickness and the local azimuthal velocity of the disk. In that case We are now able to make a boundary-layer approximation by eliminating terms involving inverse powers of the Reynolds number since here we assume that Re ) 1. Inside the boundary-layer we assume a solution to (8a)-(8h) given that the velocity components, pressure and viscosity have the following asymptotic expansions vðr;zÞ wðr;zÞ ¼ w 0 ðr; zÞ þ Re À2=ðqþ1Þ w 1 ðr; zÞ þ Á Á Á ; ð10cÞ lðr;zÞ ¼l 0 ðr; zÞ þ Re À2=ðqþ1Þl wherez ¼ z Ã =l Ã ¼ Re À1=ðqþ1Þ z is the coordinate corresponding to the region outside of the boundary-layer. To leading order we see that p 0 ¼ p 0 ðrÞ, however, since it is assumed that inside the boundarylayer the pressure is a function of z only we have that p 0 0. Having isolated the dominant viscous terms via the boundarylayer approximation, we arrive at the leading order equations that must be solved in order to determine the steady mean flow relative to the disk 1 r where the viscosity functions are given by To solve for the basic flow inside the boundary-layer we introduce the modified von Kármán [1] similarity solution in the form where g ¼ r ð1ÀqÞ=ðqþ1Þ z, thus for Bingham plastic and Carreau fluids g reduces to z. Substituting this form for the similarity solution into (11) we find that u; v; w and p must satisfy where the primes denote differentiation with respect to g and Power-law model À Bingham model À Carreau model À The model dependent form of m is given in the Appendix B along with results for the pressure correction term, p, which, as noted by Denier and Hewitt [8], is determined a posteriori once the velocity components have been calculated. By rearranging (13b) and (13c) we isolate terms involving u 00 and v 00 ; having such expressions allows us to formulate a system of five first order ordinary differential equations in five unknowns u; v; Solving (13) subject to (14) requires a full numerical solution. We employ a shooting method that utilises a fourth-order Runge-Kutta quadrature routine to perform the numerical integration of the differential equations twinned with a Newton iteration scheme to determine the values of the unknowns uð0Þ ¼ u 0 and vð0Þ ¼ v 0 .

Steady mean flow
In much the same way as Denier and Hewitt [8], we firstly consider the asymptotic form of the functions u; v and w as g ! 1.
In the limit as g ! 1 we have that w 1 For shear-thinning power-law fluids Denier and Hewitt [8] determined that ð u; vÞ $ g n=ðnÀ1Þ as g ! 1 and that n=ðn À 1Þ < 1 to ensure w is bounded in the far-field. Thus, in order to produce boundary-layer solutions that match to an external flow we are only able to consider shear-thinning fluids with n in the range 0:5 < n < 1. For shear-thickening fluids the analysis is somewhat more involved. Denier and Hewitt [8] showed that the solutions for n > 1 become non-differentiable at a critical location g ¼ g c .
However, it is noted that in this case the singularity can be completely regularised and thus the solutions can be matched to an external flow. We will not consider shear-thickening power-law fluids here, the interested reader is referred to Denier and Hewitt [8] for the full details regarding flows with n > 1 and for a lengthy discussion on shear-thinning flows with 0 < n < 0:5; n ¼ 0:5 and 0:5 < n < 1.
For Bingham plastic fluids we find that ð u; vÞ $ expð w 1 gÞ as g ! 1, as is the case for Newtonian fluids, see Cochran [11]. This is particularly interesting to note as it indicates that the large g form of the solutions has no specific dependence on the Bingham number whatsoever. The most important consequence of this result is that the velocity solutions will remain bounded as g ! 1 for all values of the Bingham number, provided, of course, that w 1 6 0. Utilising the fact that ð u 0 ; v 0 Þ ( 1 as g ! 1 we find that ð u; vÞ $ exp½ w 1 g=ð1 þ c 0 Þ for Carreau fluids as g ! 1. We note here that the exponential decay of the velocity functions has no dependence on either k or the power-law index and, more importantly, the solutions are bounded for all shear-thinning values of n, whilst for n > 1 no such 'finite width' crisis is encountered whereby the solutions become non-differentiable at a critical location, as was first noted by Dabrowski [12]. We reiterate that these solutions will exhibit exponential decay only when w 1 6 0. However, for the cases considered within this study, with c 0 ¼ 1 and k ¼ 100; 300; 500, we find that w 1 is in fact always negative for 0 < n < 1, whilst w 1 is monotonically decreasing as n increases above unity. Hence these solutions are indeed always bounded as g ! 1. Having determined the large g form of the velocity functions we are then able to apply suitable asymptotic conditions as g ! 1 within the numerical integration scheme, ensuring the correct far-field decay. This is achieved via the introduction of far-field mixed boundary conditions that relate u and v to their respective first derivatives. Solutions of (13) subject to (14) for a range of values of n; B n and k are presented in Tables 1 and 2 and Figs. 1-4, where the value of c 0 for Carreau fluids is held fixed at unity, thus making this study consistent with that of Dabrowski [12]. Where not stated, the value of 'infinity' is taken to be g 1 ¼ 20 in all cases, excluding the power-law analysis, this is sufficient to provide suitably converged results where the aforementioned asymptotic conditions are satisfied to within the desired tolerance of Oð10 À10 Þ.

Discussion
The solutions presented in Fig. 1 are identical to those of Denier and Hewitt [8] and as such we forgo any detailed analysis, noting only that the value of g 1 employed within the numerical scheme for n ¼ 0:7 and n ¼ 0:6 produces differing values of À wðg 1 Þ, for the respective values of n, when compared to the results of Denier and Hewitt [8]. This is as to be expected and is due to the ever slowing decay of the velocity functions as n approaches the critical limit of n ¼ 0:5. These results are included in this study for completeness whilst also serving as an aid to compare shearthinning results owing from the Carreau model. Plots of the viscosity function, l, are also included in Fig. 1; the growth of these functions appears to be linear with respect to the similarity coordinate. This observation is easily derived from the asymptotic analysis contained within the previous section. Given that u ¼ Ag n=ðnÀ1Þ þ Á Á Á and v ¼ Bg n=ðnÀ1Þ þ Á Á Á as g ! 1, where A and B are constants, we find that The results for Bingham plastic fluids are presented in Fig. 2. We observe a significant reduction in the peak of the cross-flow velocity component, u, as B n increases from zero, whilst the component of the azimuthal velocity, v, increases in absolute value with the Bingham number. We find that the von Kármán pumping rate, À w 1 , is decreased for increasing values of the yield stress, this being a direct consequence of the reduction in the peak of the radial velocity profile. Since, in this case, the velocity functions decay to the far-field exponentially one finds that Table 1 Numerical values of the basic flow parameters for Newtonian, power-law, Bingham plastic and Carreau fluids. The value of g 1 employed for each power-law calculation is included, we note that these are not necessarily the same as used by Denier and Hewitt [8]. Solutions for n ¼ 1 are also included for Carreau fluids, these differ from the Newtonian results as in this case the viscosity function is effectively set to l ¼ 1 þ c0, rather than unity.
Newtonian fluids u 0 where A and B are equivalent constants. Therefore for B n > 0 the viscosity functions grow exponentially in the far-field, given that w 1 6 0, as observed in Fig. 2. This unphysical result owes from the inability of the Bingham model to describe apparent viscosity at vanishing shear-rates, as noted by Zhu et al. [13].

Comparative results
A number of the previous comments regarding Bingham plastic fluids have also been noted by AS who considered the full system of non-linear governing equations by numerically integrating a three parameter system. Because of the boundary-layer formulation of this problem the governing Eq. (13) introduced here are reduced to a one parameter system, dependent only on B n . In   Re and r, thus AS must specify values for all of the parameters before a numerical solution can be obtained.
However, recalling (9a) we see that B n can be expressed simply as a function of r; Re and B y . Hence we are able to construct comparative solutions given the data used by AS, that being r ¼ 1 v; w and l versus g for Carreau fluids with n ¼ 1; 1:05; 1:25; 1:5; 1:75 and k ¼ 100. 1 Although the notation used here is not consistent with that of AS the results are directly comparable.
and Re ¼ 2950. AS have obtained solutions for B y ¼ 0; 10; 20; 30; 50 (for r ¼ 1 and Re ¼ 2950) in a stationary frame of reference. In order to produce comparative results we transform our system from a rotating reference frame to a stationary one and solve the resulting equations.
Here we have reproduced the numerical solutions of AS using their three parameter numerical scheme. The results of these calculations are presented in Fig. 5. The solid line curves are the solutions obtained from the work of AS whilst the Ã markers represent the solutions owing from the boundary-layer formulation of the problem. We note here that both sets of solutions decay to the far-field exponentially, as was outlined previously for the one parameter system. It is clear from Fig. 5 that our results are in excellent agreement with those of AS. We find that there is in fact negligible discrepancy between the two sets of solutions. This serves to confirm the boundary-layer approximation made previously. We see that the elimination of the higher order viscous terms from the problem has little to no effect on the steady mean flow solutions.
In Figs. 3 and 4 basic flow solutions are plotted for shear-thinning and shear-thickening Carreau fluids, respectively. For both cases we find that the solutions have fully converged within the confines of the boundary-layer region 0 6 g 20, which is in stark contrast to the results for shear-thinning power-law fluids presented in Fig. 1. As the power-law index increases from n ¼ 0:25 we observe that the peak in the radial cross-flow profile is shifted along the g-axis, indicating that the boundary-layer thickness increases with n. As such, the von Kármán pumping rate also increases with increasing n and does so in a non-linear fashion, see Fig. 6. We note the quantitative difference between powerlaw and Carreau shear-thinning solutions.
Clearly as g ! 1 we have that l ! 1 þ c 0 for all values of n and k. Thus, the Carreau model predicts a finite value of viscosity for both shear-thinning and shear-thickening fluids, where, in this case, the limit 1 þ c 0 is approached from below and above, respectively. Furthermore, unlike for power-law and Bingham plastic fluids, we are able to accurately model the variation of the viscosity throughout the entirety of the boundary-layer. As noted by Dabrowski [12] this is due to the unrestricted ability of the Carreau model to describe non-Newtonian fluid behaviour for both _ c Ã ( 1 and _ c Ã ) 1. Fig. 7 shows the effect of the 'relaxation' parameter, k, on the axial velocity and viscosity profiles. For shear-thinning (n ¼ 0:5) fluids, increasing the value of k has a rather minimal effect on the steady mean flow, this is observed in Table 2 as well as Fig. 7. Conversely, we see that increasing the value of k has rather a significant effect on both the velocity and viscosity profiles for shear-thickening (n ¼ 1:5) fluids. In this case both the von Kármán pumping rate and the boundary-layer thickness increases with k. Similar qualitative observations are noted by Dabrowski [12] in his Falkner-Skan boundary-layer study.

Conclusion
Solutions for the boundary-layer flow on a rotating disk have been investigated for a number of generalised Newtonian fluid models. By firstly applying a large Reynolds number boundarylayer approximation, followed by the introduction of a von Kármán-like similarity solution, the governing partial differential equations are reduced to a set of five first-order ordinary differential equations. These are solved using a fourth order Runge-Kutta quadrature routine twinned with a Newton iteration scheme. The results for shear-thinning power-law fluids are consistent with those of Denier and Hewitt [8], whilst for Bingham plastic fluids we find that implementation of the boundary-layer approximation has a negligible effect on the solutions when compared to the full-field results of AS.
Interestingly, we note that the base flow similarity solution is still applicable when considering fluids with a more complex constitutive viscosity relationship, such as the Carreau model. In this case the structure of the flow is governed by the viscosity ratio c 0 , the 'relaxation' parameter k and the power-law index n. As noted in Section 4 the Carreau model does not suffer from any of the non-physical flaws one encounters when employing either the power-law or Bingham plastic viscosity models, and as such is applicable at all shear-rates. Therefore, we suggest that the results presented within Section 3 provide a much better representation of the nature of the boundary-layer flow for both shear-thickening and shear-thinning non-Newtonian fluids. Dabrowski [12] reached a very similar conclusion in the context of flat plate boundary-layer flows.
Griffiths et al. [14,15] have indicated that shear-thinning fluids may have a stabilising effect on the rotating disk boundary-layer flow. A natural extension of this work would be to see if such results are reproduced when considering the Carreau fluid model as opposed to the simpler power-law model. Has this predicted stabilising nature been simply a consequence of the failings of the power-law model at low shear-rates?
The stability of the flow for Bingham plastic fluids could also be investigated, at least asymptotically, as in the study of Griffiths et al. [14]. However, any attempt to implement the numerical integration scheme outlined in Griffiths et al. [15] would prove futile as this scheme integrates from the outer-edge of the boundary-layer down towards the surface of the disk. Here we have shown that the viscosity profiles grow exponentially in the farfield and thus the aforementioned numerical scheme would fail to initiate. In this case at least, a differing methodology would have to be adopted in order to construct curves of neutral stability. Indeed initial investigations into both of these problems has begun; reports on these studies are to be expected in due course.
It must be stated that in the absence of experimental validation the results presented in this study must be considered as theoretical predictions only. To the best of the author's knowledge no such experiments have yet taken place, suggesting that this is an area that requires future investigation.
The one-dimensional theoretical coating investigation of Jenekhe and Schuldt [16] may be of some interest here. In this study results are presented for film thickness profiles on a rotating disk for Newtonian, power-law and Carreau fluids. The authors conclude that the breakdown of the power-law model at vanishing shear-rates has a significant effect on the predicted film thickness. As with this study the results for power-law and Carreau fluids are in stark contrast to each other. This further suggests that the applicability of the power-law model must be questioned in the context of rotating flows of this nature. As noted in Section 2 the pressure correction term, p, is realisable once the velocity functions have been calculated. The profiles in Fig. B.8 have been determined using a standard trapezoidal numerical integration scheme where the unit step-size has been reduced until sufficiently converged results were obtained. In the Newtonian limit this scheme was validated against the analytic solution, p À p 0 ¼ À2 u À w 2 =2, here we stipulate that pðg ¼ 0Þ ¼ 0 ) p 0 ¼ 0. Results for power-law and Carreau fluids are presented in Table B.3 and Fig. B.8. Due to the exponential growth of the viscosity function for Bingham plastic fluids we find that p is unbounded within the confines of the boundary-layer for B n > 0, as observed in   Numerical values of the pressure function at the outer-edge of the boundary-layer for Newtonian, power-law and Carreau fluids. The value of g 1 employed for each powerlaw calculation is included. Solutions for n ¼ 1 are also included for Carreau fluids, these differ from the Newtonian results as in this case the viscosity function is effectively set to l ¼ 1 þ c0 ¼ 2, rather than unity. Where not otherwise stated g 1 ¼ 20.

Newtonian fluids
Power-law fluids: n ¼