- Email: [email protected]
The Dean instability for shear-thinning fluids
Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
Contents lists available at SciVerse ScienceDirect
Journal of Non-Newtonian Fluid Mechanics journal homepage: http://www.elsevier.com/locate/jnnfm
The Dean instability for shear-thinning fluids Philip E. Haines a,1, James P. Denier b,⇑, Andrew P. Bassom c a
School of Mathematical Sciences, The University of Adelaide, Australia Department of Engineering Science, The University of Auckland, Auckland, New Zealand c School of Mathematics & Statistics, The University of Western Australia, Crawley 6009, Australia b
a r t i c l e
i n f o
Article history: Received 16 January 2013 Received in revised form 8 May 2013 Accepted 15 May 2013 Available online 28 May 2013 Keywords: Dean vortices Shear-thinning fluid Bifurcation analysis
a b s t r a c t We investigate the Dean instability for a generalised Newtonian fluid which satisfies an approximately power-law viscosity model, albeit modified to incorporate a low-shear Newtonian plateau. Infinite aspect ratio linear stability results are presented for both a narrow-gap width and a finite radius of curvature. These results reveal a surprising sensitivity to the details of the low-shear Newtonian region. Finite element solutions of the axisymmetric Navier–Stokes equations for flow through a finite aspect ratio duct confirm this sensitivity and, in addition, demonstrate the potential for hysteresis on the primary branch of vortices. A detailed bifurcation analysis over a range of the aspect ratio reveals that the nonlinear structure of the problem is qualitatively similar to that for a Newtonian fluid despite the apparently quite distinctive behaviour when a comparison is made at a fixed aspect ratio. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Non-Newtonian viscosity models have been successfully applied to a description of Taylor–Couette flow. In particular, for a generalised Newtonian fluid Miller and Goddard [1] demonstrated that the well-known shear-thinning index n has a significant effect on the stability characteristics of a narrow-gap Taylor flow beyond the mere adjustment to the reference viscosity. This behaviour is perhaps surprising given that, under the narrow-gap width assumption, circular Couette flow possesses a constant shear rate. The ameliorating effect can be traced to the presence of an additional multiplier of one of the radial derivatives within the governing system of equations. This extra term represents the so-called tangent viscosity, defined as the rate of change of shear stress with respect to strain rate in the direction of the base flow shear. This modification led Miller and Goddard [1] to conclude, at least for a narrow-gap width and inner cylinder rotation, that both the critical Taylor number and the critical wavenumber decrease monotonically with the index n. In contrast to the Taylor–Couette apparatus, the shear rate within curved duct flow varies significantly across the channel and vanishes at some point in the interior of the domain. One might anticipate that shear thinning has a more significant effect on the stability of this type of flow. Moreover the zero-shear phe⇑ Corresponding author. Tel.: +64 99237910. E-mail addresses: [email protected] (P.E. Haines), [email protected] (J.P. Denier), [email protected] (A.P. Bassom). 1 Present address: Department of Mathematics, University College London, United Kingdom. 0377-0257/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jnnfm.2013.05.004
nomenon provides an opportunity to study the importance of the assumptions underpinning the viscosity model at very low shear rates. We remark that the power-law model for a shear-thinning fluid (n < 1) has a viscosity singularity wherever the shear vanishes. Although the Dean instability within shear-thinning fluids has received significantly less attention than in the case of Taylor–Couette flow, a non-Newtonian viscosity model has nevertheless been applied to the Dean problem on a handful of previous occasions. The first work in this area appears to be that of Shanthini and Nandakumar [2] who examined the emergence of a four-vortex flow within a square channel. These results were soon followed by the comprehensive study of Winters [3] who fully described the bifurcation structure for the flow of a Newtonian fluid within the same geometry. It should be noted that Shanthini and Nandakumar [2] found that the critical Dean number decreases as the degree of shear thinning increases. Fellouah et al. [4] considered the stability of a power-law fluid in two particular ducts, one square and the other of aspect ratio eight. They used a Reynolds number based upon the work of Delplace and Leuliet [5], with a reference viscosity related to the average shear stress at the walls of a straight duct of identical crosssection. They reported that increasing the degree of shear thinning led to an increased critical Dean number for both aspect ratios. This apparent disagreement with the results of Shanthini and Nandakumar [2] highlights the importance of the choice of reference viscosity, a subject which we will return to in Section 3.2. In a subsequent study Fellouah et al. [6] conducted experiments using a number of shear-thinning fluids and compared their results with the earlier numerical predictions. They described the difficulty of identifying
126
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
critical Dean numbers for a developing flow confined to the 180° arc of their apparatus. Despite this inconvenience they proposed a method for identifying vortex activity at a given channel station and reported satisfactory agreement in the critical Dean number between experiment and computations. In the present study we examine the Dean instability for a shear-thinning fluid, assuming fully developed flow throughout. The implication of this assumption is that only solutions to the axisymmetric Navier–Stokes equations need be considered and this allows us to focus on the effect of introducing a shear-thinning viscosity on the fundamental development of Dean vortices. Presently, in Section 2, we describe our numerical methods and follow this in Section 3 with a description of the linear stability problem. Some finite aspect ratio results are discussed in Section 4 and the large aspect ratio case is considered in Section 5. We conclude with a brief discussion together with some final remarks.
2.1. Governing equations We consider fully developed flow through a curved channel b and with a rectangular cross-section with a radius of curvature R of width a and height b so that the aspect ratio C = b/a. It is most natural to describe the flow in terms of standard polar coordinates ð^r ; h; ^zÞ centred on the axis of curvature of the channel and as sketched in Fig. 1. We introduce a new coordinate ^ x in the radial direction such that
with
a a 6 ^x 6 ; 2 2
and work with the axisymmetric form of the Navier–Stokes equations for a fluid of constant density q. It is convenient to express the viscous components in terms of the deviatoric stress tensor so that
^ @u
^ u
^ @w
srr ¼ 2l^ ^ ; shh ¼ 2l^ ^ ; szz ¼ 2l^ ^ ; @x r @z @ v^ v^ @ v^ srh ¼ shr ¼ l^ ^ ^ ¼ ^rl^ ^ ^ ; @x r ^ @w ^ @u srz ¼ szr ¼ l^ ^ þ ^ ; @z @x @ v^ shz ¼ szh ¼ l^ ^ ; @z
a l V ^ l V ^ ¼ Vu ^t ¼ t p ^¼ p sij ¼ ð^r ; ^zÞ ¼ aðr; zÞ u sij ; V a a the Navier–Stokes equations for axisymmetric flow become
@u @u @u v 2 @p 1 @ @ srz ¼ þ Re þu þw ðr srr Þ þ @t @x @z @x r @x r @z
shh r
ð1aÞ
;
@v @v @ v uv 1 @p 1 @ @ shz ¼ Re þu þw þ þ ðr shr Þ þ r @h r @x @t @x @z r @z þ
srh r
;
ð1bÞ
@w @w @w @p 1 @ @ szz ¼ þ Re þu þw ðrszr Þ þ : @t @x @z @z r @x @z
2. Numerical method
b ^x ¼ ^r R
as this facilitates the introduction of a shear-dependent viscosity model in Section 2.2. After a nondimensionalisation in terms of : the channel width a, a velocity scale V and a reference viscosity l
@x r
ð1cÞ
Here we have retained the azimuthal pressure gradient term but discarded all other azimuthal derivatives; moreover note that b r = b + x, where b ¼ R=a is the dimensionless radius of curvature, and that the Reynolds number has been defined to be
Re ¼
qaV : l
For incompressible flow the system (1a)–(1c) must be satisfied together with
1 @ @w ðruÞ þ ¼ 0; r @r @z
ð2Þ
which is simply the conservation of mass equation for an axisymmetric flow. 2.2. The viscosity model One of the most common models used to describe fluids with a shear-dependent viscosity is the ubiquitous power-law model:
l^ ðc^Þ ¼ qjc^n1 ; ^ is the second invariant of the rate-of-strain tensor where c
c^2 ¼ 2^eij ^eij : For the axisymmetric flows considered below we have:
2 2 2 2 ^ ^ ^ ^ @w ^ 2 @u @w @ v^ v^ @u u þ2 þ2 þ þ þ ^r @ ^x ^r @ ^x @ ^z @^z @ ^x 2 @ v^ þ : @ ^z
c^2 ¼ 2
The power-law model is often applied to simple shear flows of the ^ ¼ ðuðyÞ; 0; 0Þ for which there is a single non-zero contribuform u tion to the rate of strain
@ u ^ c^ ¼ ^ : @y
Fig. 1. The curved channel geometry showing a cross-section of a channel of width a in the radial (^r or ^ x) direction and height b in the axial ð^zÞ direction. The radius of b and the aspect ratio C = b/a. curvature is R
For many such flows, for example plane Poiseuille flow, a numerical implementation of the power-law viscosity for shear-thinning fluids (0 < n < 1) must address the viscosity singularity present where ^ ¼ 0. The infinite aspect this single contribution vanishes and so c ratio Dean problem, for which the azimuthal velocity is a function of radius alone, is precisely of this type. In practice viscosity singularities tend to be mathematical artefacts of the precise form of viscosity law adopted rather than a true physical feature and in what follows we circumvent the problems presented by singularites by
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
(a)
127
(b)
Fig. 2. Effective viscosity against shear rate for both the modified power-law (5) and the Carreau (6) models and for n = 0.75 (solid line) 0.5 (dashed line) and 0.25 (dotted line). (a) Modified power-law viscosities with d = 0.01. (b) Carreau viscosities with k = 10.
considering a modified viscosity profile. Thus we adapt the powerlaw viscosity dependence and replace it with a relationship of the type n1
l^ ¼ qj½c^2 þ ^d 2 ;
ð3Þ
d is assumed to be small and clearly where the additional parameter ^ acts to smooth out the singularity at zero shear. This viscosity formulation shares features with the Carreau model [7] given by 2
n1 2
l^ ¼ l^ 1 þ ðl^ 0 l^ 1 Þ½1 þ ð^kc^Þ ; ^ 1 =l ^ 0 is often of the order 103 or 104 thereby for which the ratio l ^ 1 to be neglected in numerical calculations (e.g. [8]). This allowing l simplification then reduces the Carreau model to a three-parameter structure which is equivalent to (3) for
qj ¼ l^ 0 ^kn1 and ^d ¼ ^k2 :
ð4Þ
Somewhat inconveniently, the power-law and Carreau models tend to be nondimensionalised using different reference viscosities. In the power-law, or modified power-law, case the reference viscosity is
n1 V ; a
l^ 0 which may (but often will not) actually be realised somewhere
l ¼ qj
which is the viscosity of a power-law fluid when the dimensionless ^=V, is equal to one. This leads to a dimensionrate-of-strain, c ¼ ac less viscosity n1
l ¼ ½c2 þ d 2 ; d ¼
a2 ^ d: V2
ð5Þ
Fig. 2a depicts the effective viscosity against shear rate for a modified power-law fluid with d = 102 and several values of n < 1. On the other hand, the Carreau model is often nondimensiona¼l ^ 0 . When l ^ 1 is nelised in terms of the zero shear viscosity l glected what remains is n1
l ¼ ½1 þ ðkcÞ2 2
Fig. 3. Comparison of the modified power-law and Carreau models, with n = 0.5 in both cases, showing the respective effects of varying d and k. Results are for d = 104 (+), 102 () and 1 (solid line) and for k = 1 (solid line), 10 (dashed line) and 100 (dotted line).
with k ¼
V^ k: a
ð6Þ
Hereafter we shall only have need to refer to this reduced version of the Carreau model and so will use the term ‘Carreau’ somewhat loosely and be understood to mean a viscosity of the form (6). Fig. 2b shows the effective viscosity against shear rate for a Carreau fluid with k = 10 and a range of n < 1. We argue that the use of our proposed viscosity model has certain advantages over the standard Carreau form. Perhaps the strongest motivation for this approach is that our eventual aim is to obtain results that can be easily related to the widely used power-law case. In addition, we note that the Reynolds number is defined in terms of the effective viscosity for an order one dimensionless shear, as opposed to the theoretical maximum value
in the flow. Last, we point out that Fig. 3 demonstrates the effect of varying d within the modified power-law model and of varying k for the Carreau model. It can be seen that the influence of d is effectively confined to regions of low shear whereas k affects the viscosity at all shear rates. The conclusion is that for order one dimensionless shear rate the modified power-law viscosity is essentially determined by the single parameter n whereas the Carreau viscosity is sensitive to the parameter k as well. We remark that, to the extent that these viscosity models provide an accurate representation of fluid rheology, ^ d and ^ k should be assumed fixed for a given fluid. Thus specifying either d or k in a numerical study risks being at odds with a laboratory realisation since ^ d and ^ k are then dependent on V and so will vary with the Reynolds number. However since d influences the viscosity only over a relatively limited range of shear rates we suggest that fixing d (rather than ^ d) in a numerical investigation should be of less concern than fixing k. Finally, we mention one potential disadvantage of adopting the modified power-law approach. Its use does require the introduction of the parameter j with non-physical units and complicates the definition of the Reynolds number since
Re ¼
V 2n an
j
:
ð7Þ
Nevertheless, as the widely-used power-law model has exactly the same drawbacks, we do not consider this to be a major disadvantage here.
128
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
3. Linear theory In order to validate our finite-element computations against results for an infinite aspect ratio channel we first revisit the wellknown system of ordinary differential equations describing the linear stability of a normal mode disturbance to axisymmetric flow driven by an azimuthal pressure gradient (see, for example, Drazin ðxÞ subject to a neutral perturbation and Reid [9]). For a basic flow v of the form
Þ þ eðu ~ ; v~ ; w; ~ p ~ÞE; ðu; v ; w; pÞ ¼ ð0; v ; 0; p ~ , etc. are functions of x only, the Navier–Stokes where E eiaz and u equations can be reduced to
a2 ðDDþ þ a2 Þ½l0 ðDDþ þ a2 ÞU 2Dðl0 Dþ UÞ 2Dþ ðl0 DUÞ 2l0 ~ u bv V; ¼ 2Dn2 r þ
Fig. 5. Neutral curves for the narrow gap Dean problem with d = 102 and n as shown.
r
Dþ ðlt D VÞ þ
lt r
ð8aÞ 2
D V a
l0 V ¼ ðDþ v ÞU;
ð8bÞ
~ ; v~ Þ ¼ ðU; ReVÞ. Here we have introduced the differential where ðu operator notation
D¼
d ; dx
Dþ ¼
d 1 d 1 þ ; D ¼ ; dx r dx r
together with the nondimensional parameter, the Dean number, defined as
Dn2 ¼ Re2 b1 : Within Eq. (8b) the quantities l0 and lt are, respectively, the viscosity associated with the basic flow and the tangent viscosity: n1
n3
l0 ¼ ð½D v 2 þ dÞ 2 ; lt ¼ ðn½D v 2 þ dÞð½D v 2 þ dÞ 2 : Note that if the basic flow viscosity were constant (as is the case for circular Couette flow in the narrow-gap limit) the viscosity terms could be taken outside the derivatives. For d = 0 (so that a true power-law model is recovered) and in the limit b ? 1 the first lt term in (8b) becomes the source of the additional factor of n found by Miller and Goddard [1]. Eqs. (8a) and (8b) together with boundary conditions
1 U ¼ DU ¼ V ¼ 0 on x ¼ ; 2
ð9Þ
constitute the linear stability problem for Dean vortices in a channel with a finite radius of curvature. The narrow gap linear stability problem is recovered by taking the limit b ? 1 such that D, D+ ? D and discarding all terms involving inverse powers of r in
(8a) and (8b) excepting the ratio b/r on the right hand side of (8a) which is set to one. The system (8a), (8b) and (9) is easily solved both for finite radii of curvature and in the narrow-gap limit. We do so using a shooting method based on a fourth order Runge–Kutta scheme and, for a fixed wavenumber a, obtain a critical Dean number, Dnc, for which the boundary conditions (9) can be satisfied. On varying a we obtain the function Dnc(a) describing the location of neutral stability. 3.1. Basic flow profile Before outlining the linear stability results we discuss briefly which, in the narrowthe associated non-Newtonian base flow v gap limit, is identical to the velocity profile of an equivalent plane Poiseuille flow. Fig. 4a shows the basic radial flow profile for several values of n. As the degree of shear thinning increases so the flow profile in the central portion of the channel flattens; moreover varying the value of d does not lead to a significant change in the solution. By way of illustration for all values of n in Fig. 4a the flow profile is shown for both d = 102 and 104. In each case the two are found to overlap to within the line width used. Fig. 4b shows the effective viscosity for the same values of n. As n decreases the central portion of the flow which is at very low shear, and so large effective viscosity, widens. Altering the value of d does, however, have an appreciable effect on the viscosity in this region since it determines its maximum value. 3.2. The narrow-gap case Fig. 5 shows the neutral stability curves obtained in the limit of narrow gap width, b ? 1, for d = 102 and selected values of n. For
for d = 102 and 104. (b) Effective viscosity l for d = 102. Fig. 4. Narrow gap baseflow for n = 0.2, 0.4, 0.6, 0.8 and 1. (a) Azimuthal velocity v
129
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
each value of n 6 1 the associated neutral curve has a minimum critical Dean number, Dnc(ac) at a critical wavenumber ac(n) in the range 3 < ac < 4. As n is decreased from unity this minimum Dnc(ac(n)) initially also decreases until n is approximately 0.5 at which point Dnc starts to increase again. This behaviour is in contrast to results for the Taylor problem for which both Dnc and ac are found to decrease monotonically with the shear-thinning index n (see Miller and Goddard [1]). However, we note that in the current formulation decreasing n does not lead to an everywhere reduced viscosity as we have nondimensionalised on the effective viscosity for c = 1. In fact decreasing n lowers the viscosity where c > 1 and increases it where c < 1. Since regions of both type are present in all flows considered here it is not enough to claim that increasing the intensity of shear thinning leads to a universally lower viscosity. Fig. 6 illustrates the effect of varying d on the stability characteristics of the flow. For values of n close to unity the effect of d is quite small. Once n 6 0.6 however the situation changes somewhat and the effect of varying d is much more significant, as can be seen in Fig. 6b where varying d over two orders of magnitude, which for n = 0.4 leads to a change in the peak viscosity by a factor of four, alters Dnc(ac) by well over 10%. This result is both surprising, given that d has a negligible impact on the basic flow, and inconvenient, in that it suggests that the stability properties of shear-thinning fluids are sensitively dependent on the precise constitutive relations that they satisfy and that a single governing index is not enough to capture this behaviour. It also raises the possibility that two fluids exhibiting apparently identical flow rates and rheology may have noticeably different stability characteristics due to their behaviour at shear rates too small to be measured without specialist equipment. Given the similarity of the basic flow it might be expected that the current findings will also apply to the stability characteristics
(a)
of plane Poiseuille flow. Nouar et al. [8] and Nouar and Frigaard [10] both examined this problem for a fluid obeying the Carreau model. However, these studies only discussed the effect of k upon the strength of the shear-thinning and made no attempt to separate the effect of varying k on the low-shear Newtonian plateau. Nouar et al. [8] did, however, highlight the significance of the Reynolds number definition in interpreting the effect on stability of shear thinning. On expressing their Reynolds number in the traditional manner, that is based upon the maximum attainable viscos^ 0 , they found that increased shear thinning (decreased n or ity l increased k) led to a less stable fluid. In contrast, choosing the tangent viscosity at the wall (the minimum viscosity present anywhere in the flow) as the reference viscosity implied that shear thinning actually stabilised the flow. It should be noted that this second definition leads to a Reynolds number that depends on k. This result indicates that care must be taken when drawing conclusions regarding the effect of shear thinning upon the stability of a flow and the associated physical consequences. It is apparent from Fig. 5 that the neutral curves for n = 0.2 and n = 0.4 contain inflection points not present for larger values of n. Fig. 6b shows that this effect becomes more pronounced as d decreases. In order to investigate this change in behaviour Fig. 7 shows the neutral eigenfunctions varying with a for the two values n = 0.6 and n = 0.2. In both cases there appears to be a transition with increasing wavenumber from an eigenfunction that fills the channel to one that is confined to x > 0. This is in line with the results of Hall [11] who demonstrated that for large wavenumbers the eigenfunction is located in an asymptotic region around a critical point xc lying in the outer half of the channel. A detailed discussion of the limit a ? 1 follows in Section 3.3. For n = 0.6 the eigenfunctions appear to vary relatively smoothly as a is increased. For n = 0.2 the flow is characterised by a much larger region of very high viscosity in the centre of
(b)
65 60 55
Dn
50 45 40 35 30 2
3
4
5
α
6
7
8
9
10
2
3
4
5
α
6
7
8
9
10
Fig. 6. Effect of varying d on the neutral curves for the narrow gap width Dean problem. Results are for d = 102 (solid line), d = 103 (dashed line) and d = 104 (dotted line). (a) n = 0.8. (b) n = 0.4.
Fig. 7. Variation with a of the neutrally stable eigenfunction for the narrow gap Dean problem. (a) n = 0.6, d = 104. (b) n = 0.2, d = 104.
130
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
the channel and the transition from an eigenfunction that spans this region to one that is contained in the outer portion of the channel appears to be much more abrupt. Decreasing n increases the size of this high viscosity region whilst decreasing either n or d increases the maximum viscosity there. Both actions cause the inflection point in the neutral curve to become more pronounced and it appears that it is the rapid transition of the eigenfunction across the high viscosity central portion of the channel that is largely responsible for the unusual shape of the neutral curves reported here.
l0 ¼ jV 1 jn1 ; l1 ¼ 2l0 ðn 1ÞV 1 1 V 2; 2 2 l2 ¼ 2l0 ðn 1Þðn 2ÞV 1 V 2 þ 3l0 ðn 1ÞV 1 1 V 3: We take expansions for the perturbation as follows: 1
1
V ¼ e2 ðV 0 ðgÞ þ e2 V 1 ðgÞ þ Þ;
U ¼ U 0 ðgÞ þ e2 U 1 ðgÞ þ ;
and for the governing parameter:
1 2Dn2 ¼ e4 d0 þ e2 d1 þ :
ð13Þ
Substituting these expansions into (11) we obtain, 3.3. Large wavenumber asymptotics
1
d0 ¼ l20 ðV 0 V 1 Þ ;
Hall [11] used a large wavenumber asymptotic solution to obtain an approximation to the right hand branch of the neutral curve for both the Dean and the Görtler problems. For the Dean problem with critical parameter equal to k0 at leading order he established that the WKBJ method was valid only in regions such that 1 k0 > ðv v 0 Þ ;
ð10Þ
where a prime indicates differentiation with respect to x. When v 0 Þ1 k0 ¼ ðv min the WKBJ roots coalesce at a second-order turning point located at a critical point xc and the method breaks down. Hall [11] noted that xc is the location at which the flow violates Ray 0 j is a maxileigh’s instability criterion [9] most strongly since jv v mum there. He demonstrated that in this case solutions to the disturbance equations can still be constructed in terms of parabolic-cylinder functions in a region of thickness e1/2 about the critical point
for e = a1 a small parameter. Furthermore these solutions represent the most dangerous modes available. We now proceed to apply the method of Hall [11] to the nonNewtonian Dean problem. For simplicity we assume a narrow gap approximation and a pure power-law fluid (d = 0), since the region in which we consider expansions will be removed from the low shear plateau. In this regime the linear stability equations (8a) and (8b) reduce to:
ð15Þ
at leading order. An application of (12) at first order then yields the consistency condition d1 = 0. At second order we obtain the consistency condition:
ðn þ 2ÞD2 V 0 þ g2 3ð3 2nÞ
ð11aÞ
V1
þ ð5 2nÞ
V2
dv 2 l ; dx
and so, noting that l ¼ jv 0 jn1 , to a critical location for the instabil 0 jv 0 j22n is a minimum: ity such that v v
V0
d2 ¼ d0 ðn þ 2Þ 3ð3 2nÞ
V3 V1
þ ð5 2nÞ
V2 V0
Dn ¼ a2
1 ðd0 þ a1 d2 Þ 2
ð16Þ
:
1=2 ð17Þ
;
where d0 and d2 are obtained from the basic flow profile according to (14) and (16). 3.4. The finite-gap case Results for the properties of Dean vortices in channels with a finite radii of curvature are limited, even for Newtonian fluids. Finlay and Nandakumar [12] calculated both ac and Dnc(ac) to very high accuracy and provided some detailed comparison with the results
1000
100
n= 1
1
l ¼ l0 þ e2 gl1 þ eg2 l2 þ
0.8 0.6 0.4
where
and
!#12
Fig. 8 shows neutral curves for d = 104 and various values of n together with the two term asymptotic prediction for the neutral Dean number
v ¼ V 0 þ e gV 1 þ eg2 V 2 þ
x¼xc
d2 V 0 ¼ 0: d0
and for non-negative integer values of m. The most unstable mode is therefore:
1 2
10
;
V0 þ
!#12 " 1 V3 V2 d2 ¼ 2d0 þ m ðn þ 2Þ 3ð3 2nÞ þ ð5 2nÞ ; 2 V1 V0
ð12Þ
Note that the variable viscosity has severed the link between the critical point and the point at which Rayleigh’s criterion is most strongly violated. In order to expand the governing Eqs. (11) about the critical point we expand the basic flow and its viscosity in terms of g
j 1 d v Vj ¼ j! dxj
!
ð11bÞ
An application of the WKBJ method to these equations leads to roots
d 2 ðv v 0 jv 0 j22n Þ ¼ l2 ðð3 2nÞv v 00 þ ðv 0 Þ Þ ¼ 0: dx
V3
Solutions of this equation can be obtained in terms of the paraboliccylinder function and are such that they decay as jgj ? 1 when
Dn
e2 ðD2 þ e2 Þ½lðD2 þ e2 ÞU 4DðlDUÞ ¼ 2Dn2 v V nDðlDVÞ e2 lV ¼ ðDv ÞU: 2
U 0 ¼ l
"
1
g ¼ e2 ðx xc Þ;
ðK 2 þ 1Þ ðnK 2 þ 1Þ ¼ k0 v
ð14Þ
1 0V 1 V 0;
1
α
10
Fig. 8. Logarithmic plot of the neutral curves for various values of n and d = 104 together with the two term asymptotic prediction for the limit a ? 1 given by (17).
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
131
4. Finite aspect ratio: A numerical study 4.1. Finite element method
Fig. 9. Neutral curves varying with curvature ratio, b, for n = 0.6 and d = 102.
We now examine whether the infinite aspect ratio linear stability results of Section 3 are applicable to a domain of finite extent in the z-direction, as depicted in Fig. 1. Of particular interest is whether the stability characteristics remain sensitive to d when the base flow is no longer a function of radius alone and when disturbances of finite amplitude are admitted. The nondimensionalised governing equations for axisymmetric flow were given in Section 2.1. Here we solve (1a)–(1c) and (2) using a finite-element method. We adopt the same approach as in Haines et al. [14] (herein referred to as HDB) and so only outline the main points of the procedure. The system of equations is solved on the rectangular domain
Fig. 10. Dnc against n for a = 3.93, b = 19.5 and d = 102 (solid line) and d = 104 (dashed line).
of previous studies. Our results for n = 1 agree with those of Finlay and Nandakumar [12] to at least four decimal places in both Dnc and ac. For n – 1 our linear stability analysis confirms that, as for the Newtonian case [13], the critical Dean number is significantly greater than that for the narrow-gap approximation, even for quite large values of b. Fig. 9 shows neutral curves for n = 0.6, d = 102 and a range of curvature ratios. We suggest that only for b > 100 can the neutral curve effectively be considered independent of the radius of curvature. Critical values of the Dean number against n for a = 2p/ 1.6 = 3.93 and b = 19.5 are shown in Fig. 10. We note that Fig. 10 does not indicate the onset of instability, since, for each value of n, it gives the critical value of the Dean number for a fixed a = 3.93 rather than the minimum taken over all wavenumbers Dnc(ac(n)). The value of a chosen is intended for later use in Section 4.2 where the results for an infinite aspect ratio presented here will be compared with computations in a channel of finite aspect ratio C = 1.6. For the values of n considered the wavenumber a = 3.93 is, however, always of a similar magnitude to the critical values which lie in the range 3 < ac(n) < 4. Fig. 10 indicates that d has a significant impact on the stability of the flow for values of n away from one. It also serves to illustrate the key findings of our linear stability study which we summarise as follows:
1 1 6x6 ; 2 2
1 6 Z 6 1;
ð18Þ
where we have defined a new coordinate Z = 2z/C. This has the advantage of introducing C alongside b as an explicit parameter in the equations so that continuation in either is possible. The computational domain is sketched in Fig. 11 and the weak form used to represent 1a,1b,1c and 2 can be found in Appendix A. The computational mesh is uniformly spaced in the coordinates x and Z and comprises (initially) Nx elements in the radial direction and NZ elements in the axial direction. A typical mesh resolution was of the order Nx = 20, NZ = 40. Experiments were conducted with various mesh dimensions and it was checked that the results discussed below were obtained on a sufficiently refined grid. Moreover, if any particular calculations required finer resolution, it was possible, as necessary, to adaptively refine the mesh as part of the solution procedure. The Z2 error estimator developed by Zienkiewicz and Zhu [15] with a flux based upon the shear stress was used to control the size of the mesh elements. The finite-element library oomph-lib Heil and Hazel [16] was used to assemble and solve the system of algebraic equations associated with a weak form of the axisymmetric Navier–Stokes equations. The Q2P1 finite element (see Gresho and Sani [17], p. 554) was used throughout. 4.2. Single wavelength validation In order to validate our finite-element computations against those for an infinite aspect ratio channel we first considered the case of a channel subject to symmetry conditions
w¼
@u @ v ¼ ¼ 0; @Z @Z
Decreasing n from one initially lowers the neutral curve, and in particular decreases Dnc(ac). This trend is reversed for n [ 0.5. For n [ 0.8 the value of d significantly affects Dnc. To the best of our knowledge these findings remain valid regardless of the curvature ratio.
Fig. 11. The axisymmetric computational geometry in terms of the coordinates given in (18).
132
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
at Z = ± 1. These boundary conditions allow the computation of a single vortex pair, the emergence of which with increasing Dean number can be compared with that predicted by the infinite aspect ratio linear stability analysis. On the side walls, x = ± 1/2, all velocities are set to zero with the exception of a small azimuthal velocity (everywhere less than 105 times the peak azimuthal velocity) on the inner cylinder boundary. This serves to slightly disconnect the bifurcation at Dnc by preferentially selecting the branch of vortices with outflow along the centreline of the domain. Fig. 12 shows results for C = 1.6, b = 19.5 and selected values of n 6 1. The radial velocity at the centre of the domain u(0, 0) is used as a suitable measure of the solutions. The vertical lines in Fig. 12 indicate the position at which a loss of stability (to the specified wavenumber) is predicted by linear theory. The agreement between linear theory and the two-dimensional numerical results is found to be very good. In particular Fig. 12 confirms that the Dean number for the onset of a single vortex pair does not decrease monotonically with n. Fig. 13 reproduces Fig. 12 with additional data for d = 104. The agreement with linear theory is again found to be good, with the results in Fig. 13 confirming the increasing influence of d with decreasing n. Fig. 13 also exhibits some striking features not captured by the linear theory. Firstly, although the value of d significantly affects
Fig. 12. Solution branches for a single vortex pair of wavelength 1.6 contained in a finite domain subject to symmetry conditions at Z = ± 1. Results are for C = 1.6, b = 19.5 and d = 102. Vertical lines and accompanying values indicate the position at which a loss of stability is predicted by linear theory. The measure of the flow used is the centreline radial velocity at the point x = Z = 0.
Fig. 13. Solution branches for a single vortex pair of wavelength 1.6 contained in a finite domain subject to symmetry conditions at Z = ± 1. Results are for C = 1.6, b = 19.5 and d = 104. Results for d = 102, previously given in Fig. 12, are shown in grey.
the location of Dnc, it has less of an impact on the vortex modes once they become nonlinear. For all three values of n < 1 shown in Fig. 13 the solution branches for the two values of d move closer together as the strength of the vortices (indicated by ju(0, 0)j) increases. This can be attributed to a reduction in the amount of fluid at very low shear and to the corresponding decrease in importance of the details of the zero-shear plateau. The second striking feature of Fig. 13 is the presence of hysteresis for d = 104. This hysteresis becomes more pronounced as n is reduced. For n = 0.4 plausibly stable nonlinear vortex modes exist for Dean numbers in the range 35.2 < Dn < Dnc = 35.83. 4.3. Comparison with Shanthini and Nandakumar Shanthini and Nandakumar [2] explored the effect of the power-law index n on the bifurcation structure in a square duct. They performed symmetric computations that identified the stable two and four-cell parts of the main solution branch and noted that the non-Newtonian characteristics of the fluid did not appear to change the qualitative behaviour of the flow. The turning point at which the two-vortex solution jumps to a four-vortex solution with increasing Dn (Dn24 in Shanthini’s notation) and the corresponding point at which the four-vortex solution jumps to a two-vortex solution with decreasing Dn (Dn42) are, for n = 1, the bifurcation points L1 and L2 described by Winters [3] (see his Fig. 3). A recreation of Fig. 6 of Shanthini and Nandakumar [2] showing the location of the two bifurcation points in the Dn–n plane is given in Fig. 14. Our results, obtained by using bifurcation detection to track the turning points Dn24 and Dn42 over a range of n, are shown as solid and dashed lines respectively. The original data points of Shanthini and Nandakumar [2], computed to an accuracy of Dnc ± 1 using the bisection method, have been estimated graphically. The agreement is generally good, although there appears to be a departure for Dn24 when n = 0.6. At such small aspect ratios the flow experiences non-zero shear everywhere apart from exactly in the corners of the domain. Shanthini and Nandakumar [2] dealt with this corner singularity by imposing a maximum viscosity of 106 there, far greater than was found anywhere else in the flow. Although they do not appear to have performed any sensitivity analysis on their viscosity ceiling we have observed that for a square duct and d 6 104 the correction to a purely power-law viscosity is felt only in the corners and that any further reduction in d does not influence the bifurcation loci significantly. Thus the parameter d is found to be of
Fig. 14. Loci of fold bifurcations in the Dn–n plane for b = 100 and C = 1. Solid line: Turning point at which two vortex solution jumps to four vortex solution with increasing Dn, Dn24 in Shanthini’s notation. Dashed line: Turning point at which four vortex solution jumps to two vortex solution with decreasing Dn, Dn42 in Shanthini’s notation. Shanthini’s data points are also shown, they have been reconstructed from Fig. 6 of the original paper.
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
133
limited importance for such small aspect ratios. Consequently the power-law model can be, and has been, used without the care that we shall see is required for larger aspect ratios, for which regions of very low shear persist in the interior of the domain.
5. Large aspect ratio results We now consider the axisymmetric solution structure for curved channels of large aspect ratio, which we take to mean C 10. The corresponding structure for a Newtonian fluid was previously described in HDB. For C = 10 we solve the governing equations for flow in the upper half channel Z > 0 subject to the symmetry boundary conditions w = uZ = vZ = 0 on Z = 0. The velocities on all other boundaries are set to zero. We first compute the unique low Dean number solution and then use arc-length continuation to track the associated (primary) solution branch to higher Dean number. Fig. 15 depicts this solution branch close to the onset of instability for a modified power-law fluid with shear-thinning index n = 0.6 and d = 102. Instability is found to occur at a lower Dean number than for the Newtonian case, as suggested by the linear stability calculations in Section 3. Fig. 15 also shows streamline plots of the secondary velocity at selected points. Note that in all cases the flow is symmetric about Z = 0 and so only the half of the domain Z > 0 is shown. As for a Newtonian fluid the flow at low Dean number is characterised by a large end vortex but an absence of vortex motion on the channel centreline (u(0, 0) 0). As the Dean number is increased through the critical value for centrifugal instability the flow on the primary branch changes to one with vortices throughout the channel. Fig. 15 also illustrates two disconnected solution branches. We refer to the solutions on these branches as secondary modes. A streamline plot shows the crosssectional flow for one of these modes and it can be seen that the number of vortices in the half channel differs by one from the primary mode. The (secondary) branch of solutions characterised by outflow along the centreline is formed of two distinct parts. Whilst this was not the case for a Newtonian fluid at the same aspect ratio, similar behaviour was observed in the Newtonian case at slightly smaller aspect ratios (see HDB). We shall return to this observation in Section 5.1 below. The effect of decreasing d upon the bifurcation structure depicted in Fig. 15 is shown in Fig. 16. Reducing d by two orders of magnitude (corresponding to an increase in the fluid’s zero-shear
Fig. 16. Symmetric solution branches for n = 0.6 and for d = 102 (solid line), 103 (dashed line) and 104 (dotted line). For d = 103 and 104 we have made no attempt to continue the solution branches beyond Dn 29.4.
viscosity of roughly 2.5) leads to a change of around 4% in the critical Dean number. This is qualitatively in line with both the linear stability results of Fig. 10 and the single wavelength results of Section 4.2. We do not expect quantitative agreement for a number of reasons discussed in HDB. These include the difficulty of establishing a wavenumber for the onset of the finite aspect ratio instability and a discrepancy in the magnitude of the nondimensionalised centreline velocity between the finite and infinite aspect ratio results. Both of these phenomena may be attributed to the presence of end effects. Although an aspect ratio of 10 may be regarded as large it is clear from Fig. 15 that the end vortex region still takes up a significant portion of the channel both with and without the presence of the instability. As in Section 4.2 the results in Fig. 16 also enable the effect of d on nonlinear perturbation growth to be examined. As was noted for a single wavelength, reducing d introduces a hysteresis region to the primary solution branch. The size of this hysteresis region is observed to increase with decreasing d and also with decreasing n (although this effect is not shown here). Again, just as in Section 4.2, we find that the impact of d is at its most pronounced close to the onset of instability. As the strength of the vortices grow throughout the channel (indicated by increasing ju(0, 0)j) d has a decreasing impact on the solution branches due to a reduction in size of the low shear regions where its effect is felt.
Fig. 15. Symmetric solution branches for a modified power-law fluid with n = 0.6 and d = 102 and for C = 10, b = 19.5. The labels attached to the limit points correspond to those in Fig. 17. The inset streamline plots show the secondary flow found on selected solution branches.
134
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
Fig. 17. The location of limit points in the Dn C plane for (a) n = 0.6 and d = 102 and (b) the Newtonian fluid as reported previously in HDB. Note that the two figures cover different coordinate ranges. The labels I and II relate the sets of limit points shown in (a) to those given (for C = 10) in Fig. 15.
5.1. Bifurcation structure with varying aspect ratio Neither the solution branches or secondary streamlines depicted in Fig. 15 bear much resemblance to those for a Newtonian fluid at the same aspect ratio (see, for example, Fig. 3 of HDB). However, the rapid change of these features with aspect ratio observed by the previous study indicates that a description of the flow at a single aspect ratio does not tell the full story. Hence now we detail results for the loci of limit points in the Dn C plane. Corresponding loci for the Newtonian problem were given in HDB where it was found that the primary mode of instability (characterised by the number of vortices accommodated in the channel cross-section) was selected via a transcritical bifurcation between two different modes. This process was first described by Benjamin [18] in the context of the Taylor problem. In order to examine the non-Newtonian bifurcation structure the two turning point bifurcations shown in Fig. 15 (for n = 0.6 and d = 102) have been continued over a range of aspect ratios 6 < C < 13. The results are shown in Fig. 17a. A comparison with Fig. 6 of HDB (reproduced in part here as Fig. 17b) demonstrates that qualitatively similar features are found at an aspect ratio somewhere between one and two greater for n = 0.6 and d = 102 than in the Newtonian case. For example the bifurcation structure for C = 10 given in Fig. 15 is found to be similar to that found for a Newtonian fluid with C 8 (see the middle schematic of Fig. 7a of HDB). Although several other turning point bifurcations were found for the Newtonian problem we have made no attempt to compute them here. As Shanthini and Nandakumar [2] observed for C = 1, shear thinning does not appear to change the qualitative structure of the bifurcation set. 6. Discussion We introduced a new viscosity function (3) with similarities to both the power-law and Carreau models. In addition to the shearthinning index, n, this viscosity function contains a second parameter d that controls the location of the low-shear Newtonian plateau. This second parameter has the advantage over the Carreau model that it does not significantly affect the viscosity at dimensionless strain rates of order one. Nevertheless, and despite having a negligible impact on the basic flow profile d is found to have a significant impact upon the onset of instability. This linear stability finding is confirmed by fully nonlinear results for axisymmetric flow in a finite aspect ratio duct, computations that additionally demonstrate the existence of hysteresis on the primary branch of vortices. Decreasing either n or d increases the size of the
hysteresis region. As the nonlinear vortices increase in strength so the proportion of the domain at very low shear decreases and consequently so does the influence of d. Although the structure in Fig. 15 appears to be very different from that for a Newtonian fluid (Fig. 3 of HDB) a bifurcation analysis shows that this can be attributed to a shift in the aspect ratio at which qualitatively similar features are present rather than to a fundamental change in the character of the vortices. Given the sensitive dependence upon aspect ratio of the vortex structures present in both the Dean and Taylor–Couette problems, bifurcation tracking provides a powerful tool with which to unravel the underlying behaviour. Finally we speculate as to whether plane Poiseuille flow exhibits a similar sensitivity to the low-shear Newtonian plateau. Although the nature of the two instabilities is fundamentally different the basic flow profile in the two problems is very similar being characterised by a velocity profile that is a function of a single transverse coordinate and that has a point of vanishing shear somewhere in the interior of the flow. Whilst the instability of plane Poiseuille flow of a shear-thinning fluid has been previously studied [8,10] the use of the Carreau model meant it was not possible to isolate the effect of k upon the Newtonian plateau from its effect on the viscosity at larger shear rates. Acknowledgements JPD and APB gratefully acknowledge the support of the Australian Research Council. The comments of the anonymous referees are gratefully acknowledged. Appendix A. Weak form In the finite-element computations Eqs. (1a)–(1c) and (2) were represented by the following weak forms:
Rui ¼
Z Z
@u @u w @u v 2 þu þ Wi r dx dZ @t @x r @Z r X Z Z p u 2l 2 Wi r dx dZ r X r Z Z Z Z @u @ Wi 1 @u @w 1 p 2l l þ r dx dZ þ @x @x r @Z @x r X X @ Wi
r dx dZ; @Z St
where we define Rui ð1 6 i 6 NÞ to be the residual associated with evaluating the r-momentum equation with respect to the test function Wi, r = C/2 and X is our computational domain (18). Similarly:
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
Z Z @v @ v w @ v uv Rvi ¼ St þu þ þ Wi r dx dZ @t @x r @Z r X Z Z l @v v 1 @p Wi r dx dZ r @h r @x r X Z Z Z Z @ v v @ Wi 1 @v 1 þ l l r dx dZ þ @x r @x r @Z r X X @ Wi
r dx dZ; @Z and
Rw i ¼
Z Z @w @w w @w St þu þ Wi r dx dZ @t @x r @Z X Z Z Z Z @w 1 @u @ Wi l @w 1 l þ p2 þ r dx dZ @x r @Z @x r @Z r X X @ Wi
r dx dZ; @Z
The weak form for conservation of mass is constructed using the (linear) pressure basis functions Uj (1 6 j 6 M) as test functions:
Rpj ¼
Z Z @u u 1 @w þ þ Uj r dx dZ: r r @Z X @r
References [1] C. Miller, J. Goddard, Appendix to polymer fluid mechanics, Adv. Appl. Mech. 19 (1979) 143–219. [2] W. Shanthini, K. Nandakumar, Bifurcation phenomena of generalized Newtonian fluids in curved rectangular ducts, J. Non-Newton. Fluid Mech. 22 (1) (1986) 35–60.
135
[3] K. Winters, Bifurcation study of laminar flow in a curved tube of rectangular cross-section, J. Fluid Mech. 180 (1987) 343–369. [4] H. Fellouah, C. Castelain, A. Ould-El-Moctar, H. Peerhossaini, A numerical study of Dean instability in non-Newtonian fluids, J. Fluids Eng. 128 (2006) 34–41. [5] F. Delplace, J.C. Leuliet, Generalized Reynolds number for the flow of power law fluids in cylindrical ducts of arbitrary cross-section, Chem. Eng. J. Biochem. Eng. J. 56 (1995) 33–37. [6] H. Fellouah, C. Castelain, A. Ould-El-Moctar, H. Peerhossaini, The Dean instability in power-law and Bingham fluids in a curved rectangular duct, J. Non-Newton. Fluid Mech. 165 (3–4) (2010) 163–173. [7] R. Bird, R. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, Fluid Mechanics, vol. 1, John Wiley and Sons Inc., New York, NY, 1987. [8] C. Nouar, A. Bottaro, J. Brancher, Delaying transition to turbulence in channel flow: revisiting the stability of shear-thinning fluids, J. Fluid Mech. 592 (2007) 177–194. [9] P.G. Drazin, W.H. Reid, Hydrodynamic Stability, Cambridge University Press, 1981. [10] C. Nouar, I. Frigaard, Stability of plane Couette–Poiseuille flow of shearthinning fluid, Phys. Fluids 21 (2009) 064104. [11] P. Hall, Taylor–Görtler vortices in fully developed or boundary-layer flows: linear theory, J. Fluid Mech. 124 (1) (1982) 475–494. [12] W. Finlay, K. Nandakumar, Onset of two-dimensional cellular flow in finite curved channels of large aspect ratio, Phys. Fluids A: Fluid Dynam. 2 (1990) 1163–1174. [13] E. Sparrow, On the onset of flow instability in a curved channel of arbitrary height, Z. Angew. Math. Phys. (ZAMP) 15 (6) (1964) 638–642. [14] P.E. Haines, J.P. Denier, A.P. Bassom, The Dean instability for finite aspect ratio ducts, J. Fluid Mech. 716 (2013). R8-1-13. [15] O. Zienkiewicz, J. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique, Int. J. Numer. Methods Eng. 33 (7) (1992) 1331–1364. [16] M. Heil, A. Hazel, oomph-lib-An object-oriented multi-physics finite-element library, in: M. Schafer, H.-J. Bungartz (Eds.), Fluid-Structure Interaction, Lecture Notes on Computational Science and Engineering, 2006, pp. 19–49. [17] P. Gresho, R. Sani, Incompressible Flow and the Finite Element Method, Advection-Diffusion and Isothermal Laminar Flow, vol. 1, John Wiley & Sons, 1998. [18] T. Benjamin, Bifurcation phenomena in steady flows of a viscous fluid. I. Theory, Proc. Roy. Soc. Lond. A. Math. Phys. Sci. 359 (1696) (1978) 1–26.
Contents lists available at SciVerse ScienceDirect
Journal of Non-Newtonian Fluid Mechanics journal homepage: http://www.elsevier.com/locate/jnnfm
The Dean instability for shear-thinning fluids Philip E. Haines a,1, James P. Denier b,⇑, Andrew P. Bassom c a
School of Mathematical Sciences, The University of Adelaide, Australia Department of Engineering Science, The University of Auckland, Auckland, New Zealand c School of Mathematics & Statistics, The University of Western Australia, Crawley 6009, Australia b
a r t i c l e
i n f o
Article history: Received 16 January 2013 Received in revised form 8 May 2013 Accepted 15 May 2013 Available online 28 May 2013 Keywords: Dean vortices Shear-thinning fluid Bifurcation analysis
a b s t r a c t We investigate the Dean instability for a generalised Newtonian fluid which satisfies an approximately power-law viscosity model, albeit modified to incorporate a low-shear Newtonian plateau. Infinite aspect ratio linear stability results are presented for both a narrow-gap width and a finite radius of curvature. These results reveal a surprising sensitivity to the details of the low-shear Newtonian region. Finite element solutions of the axisymmetric Navier–Stokes equations for flow through a finite aspect ratio duct confirm this sensitivity and, in addition, demonstrate the potential for hysteresis on the primary branch of vortices. A detailed bifurcation analysis over a range of the aspect ratio reveals that the nonlinear structure of the problem is qualitatively similar to that for a Newtonian fluid despite the apparently quite distinctive behaviour when a comparison is made at a fixed aspect ratio. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Non-Newtonian viscosity models have been successfully applied to a description of Taylor–Couette flow. In particular, for a generalised Newtonian fluid Miller and Goddard [1] demonstrated that the well-known shear-thinning index n has a significant effect on the stability characteristics of a narrow-gap Taylor flow beyond the mere adjustment to the reference viscosity. This behaviour is perhaps surprising given that, under the narrow-gap width assumption, circular Couette flow possesses a constant shear rate. The ameliorating effect can be traced to the presence of an additional multiplier of one of the radial derivatives within the governing system of equations. This extra term represents the so-called tangent viscosity, defined as the rate of change of shear stress with respect to strain rate in the direction of the base flow shear. This modification led Miller and Goddard [1] to conclude, at least for a narrow-gap width and inner cylinder rotation, that both the critical Taylor number and the critical wavenumber decrease monotonically with the index n. In contrast to the Taylor–Couette apparatus, the shear rate within curved duct flow varies significantly across the channel and vanishes at some point in the interior of the domain. One might anticipate that shear thinning has a more significant effect on the stability of this type of flow. Moreover the zero-shear phe⇑ Corresponding author. Tel.: +64 99237910. E-mail addresses: [email protected] (P.E. Haines), [email protected] (J.P. Denier), [email protected] (A.P. Bassom). 1 Present address: Department of Mathematics, University College London, United Kingdom. 0377-0257/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jnnfm.2013.05.004
nomenon provides an opportunity to study the importance of the assumptions underpinning the viscosity model at very low shear rates. We remark that the power-law model for a shear-thinning fluid (n < 1) has a viscosity singularity wherever the shear vanishes. Although the Dean instability within shear-thinning fluids has received significantly less attention than in the case of Taylor–Couette flow, a non-Newtonian viscosity model has nevertheless been applied to the Dean problem on a handful of previous occasions. The first work in this area appears to be that of Shanthini and Nandakumar [2] who examined the emergence of a four-vortex flow within a square channel. These results were soon followed by the comprehensive study of Winters [3] who fully described the bifurcation structure for the flow of a Newtonian fluid within the same geometry. It should be noted that Shanthini and Nandakumar [2] found that the critical Dean number decreases as the degree of shear thinning increases. Fellouah et al. [4] considered the stability of a power-law fluid in two particular ducts, one square and the other of aspect ratio eight. They used a Reynolds number based upon the work of Delplace and Leuliet [5], with a reference viscosity related to the average shear stress at the walls of a straight duct of identical crosssection. They reported that increasing the degree of shear thinning led to an increased critical Dean number for both aspect ratios. This apparent disagreement with the results of Shanthini and Nandakumar [2] highlights the importance of the choice of reference viscosity, a subject which we will return to in Section 3.2. In a subsequent study Fellouah et al. [6] conducted experiments using a number of shear-thinning fluids and compared their results with the earlier numerical predictions. They described the difficulty of identifying
126
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
critical Dean numbers for a developing flow confined to the 180° arc of their apparatus. Despite this inconvenience they proposed a method for identifying vortex activity at a given channel station and reported satisfactory agreement in the critical Dean number between experiment and computations. In the present study we examine the Dean instability for a shear-thinning fluid, assuming fully developed flow throughout. The implication of this assumption is that only solutions to the axisymmetric Navier–Stokes equations need be considered and this allows us to focus on the effect of introducing a shear-thinning viscosity on the fundamental development of Dean vortices. Presently, in Section 2, we describe our numerical methods and follow this in Section 3 with a description of the linear stability problem. Some finite aspect ratio results are discussed in Section 4 and the large aspect ratio case is considered in Section 5. We conclude with a brief discussion together with some final remarks.
2.1. Governing equations We consider fully developed flow through a curved channel b and with a rectangular cross-section with a radius of curvature R of width a and height b so that the aspect ratio C = b/a. It is most natural to describe the flow in terms of standard polar coordinates ð^r ; h; ^zÞ centred on the axis of curvature of the channel and as sketched in Fig. 1. We introduce a new coordinate ^ x in the radial direction such that
with
a a 6 ^x 6 ; 2 2
and work with the axisymmetric form of the Navier–Stokes equations for a fluid of constant density q. It is convenient to express the viscous components in terms of the deviatoric stress tensor so that
^ @u
^ u
^ @w
srr ¼ 2l^ ^ ; shh ¼ 2l^ ^ ; szz ¼ 2l^ ^ ; @x r @z @ v^ v^ @ v^ srh ¼ shr ¼ l^ ^ ^ ¼ ^rl^ ^ ^ ; @x r ^ @w ^ @u srz ¼ szr ¼ l^ ^ þ ^ ; @z @x @ v^ shz ¼ szh ¼ l^ ^ ; @z
a l V ^ l V ^ ¼ Vu ^t ¼ t p ^¼ p sij ¼ ð^r ; ^zÞ ¼ aðr; zÞ u sij ; V a a the Navier–Stokes equations for axisymmetric flow become
@u @u @u v 2 @p 1 @ @ srz ¼ þ Re þu þw ðr srr Þ þ @t @x @z @x r @x r @z
shh r
ð1aÞ
;
@v @v @ v uv 1 @p 1 @ @ shz ¼ Re þu þw þ þ ðr shr Þ þ r @h r @x @t @x @z r @z þ
srh r
;
ð1bÞ
@w @w @w @p 1 @ @ szz ¼ þ Re þu þw ðrszr Þ þ : @t @x @z @z r @x @z
2. Numerical method
b ^x ¼ ^r R
as this facilitates the introduction of a shear-dependent viscosity model in Section 2.2. After a nondimensionalisation in terms of : the channel width a, a velocity scale V and a reference viscosity l
@x r
ð1cÞ
Here we have retained the azimuthal pressure gradient term but discarded all other azimuthal derivatives; moreover note that b r = b + x, where b ¼ R=a is the dimensionless radius of curvature, and that the Reynolds number has been defined to be
Re ¼
qaV : l
For incompressible flow the system (1a)–(1c) must be satisfied together with
1 @ @w ðruÞ þ ¼ 0; r @r @z
ð2Þ
which is simply the conservation of mass equation for an axisymmetric flow. 2.2. The viscosity model One of the most common models used to describe fluids with a shear-dependent viscosity is the ubiquitous power-law model:
l^ ðc^Þ ¼ qjc^n1 ; ^ is the second invariant of the rate-of-strain tensor where c
c^2 ¼ 2^eij ^eij : For the axisymmetric flows considered below we have:
2 2 2 2 ^ ^ ^ ^ @w ^ 2 @u @w @ v^ v^ @u u þ2 þ2 þ þ þ ^r @ ^x ^r @ ^x @ ^z @^z @ ^x 2 @ v^ þ : @ ^z
c^2 ¼ 2
The power-law model is often applied to simple shear flows of the ^ ¼ ðuðyÞ; 0; 0Þ for which there is a single non-zero contribuform u tion to the rate of strain
@ u ^ c^ ¼ ^ : @y
Fig. 1. The curved channel geometry showing a cross-section of a channel of width a in the radial (^r or ^ x) direction and height b in the axial ð^zÞ direction. The radius of b and the aspect ratio C = b/a. curvature is R
For many such flows, for example plane Poiseuille flow, a numerical implementation of the power-law viscosity for shear-thinning fluids (0 < n < 1) must address the viscosity singularity present where ^ ¼ 0. The infinite aspect this single contribution vanishes and so c ratio Dean problem, for which the azimuthal velocity is a function of radius alone, is precisely of this type. In practice viscosity singularities tend to be mathematical artefacts of the precise form of viscosity law adopted rather than a true physical feature and in what follows we circumvent the problems presented by singularites by
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
(a)
127
(b)
Fig. 2. Effective viscosity against shear rate for both the modified power-law (5) and the Carreau (6) models and for n = 0.75 (solid line) 0.5 (dashed line) and 0.25 (dotted line). (a) Modified power-law viscosities with d = 0.01. (b) Carreau viscosities with k = 10.
considering a modified viscosity profile. Thus we adapt the powerlaw viscosity dependence and replace it with a relationship of the type n1
l^ ¼ qj½c^2 þ ^d 2 ;
ð3Þ
d is assumed to be small and clearly where the additional parameter ^ acts to smooth out the singularity at zero shear. This viscosity formulation shares features with the Carreau model [7] given by 2
n1 2
l^ ¼ l^ 1 þ ðl^ 0 l^ 1 Þ½1 þ ð^kc^Þ ; ^ 1 =l ^ 0 is often of the order 103 or 104 thereby for which the ratio l ^ 1 to be neglected in numerical calculations (e.g. [8]). This allowing l simplification then reduces the Carreau model to a three-parameter structure which is equivalent to (3) for
qj ¼ l^ 0 ^kn1 and ^d ¼ ^k2 :
ð4Þ
Somewhat inconveniently, the power-law and Carreau models tend to be nondimensionalised using different reference viscosities. In the power-law, or modified power-law, case the reference viscosity is
n1 V ; a
l^ 0 which may (but often will not) actually be realised somewhere
l ¼ qj
which is the viscosity of a power-law fluid when the dimensionless ^=V, is equal to one. This leads to a dimensionrate-of-strain, c ¼ ac less viscosity n1
l ¼ ½c2 þ d 2 ; d ¼
a2 ^ d: V2
ð5Þ
Fig. 2a depicts the effective viscosity against shear rate for a modified power-law fluid with d = 102 and several values of n < 1. On the other hand, the Carreau model is often nondimensiona¼l ^ 0 . When l ^ 1 is nelised in terms of the zero shear viscosity l glected what remains is n1
l ¼ ½1 þ ðkcÞ2 2
Fig. 3. Comparison of the modified power-law and Carreau models, with n = 0.5 in both cases, showing the respective effects of varying d and k. Results are for d = 104 (+), 102 () and 1 (solid line) and for k = 1 (solid line), 10 (dashed line) and 100 (dotted line).
with k ¼
V^ k: a
ð6Þ
Hereafter we shall only have need to refer to this reduced version of the Carreau model and so will use the term ‘Carreau’ somewhat loosely and be understood to mean a viscosity of the form (6). Fig. 2b shows the effective viscosity against shear rate for a Carreau fluid with k = 10 and a range of n < 1. We argue that the use of our proposed viscosity model has certain advantages over the standard Carreau form. Perhaps the strongest motivation for this approach is that our eventual aim is to obtain results that can be easily related to the widely used power-law case. In addition, we note that the Reynolds number is defined in terms of the effective viscosity for an order one dimensionless shear, as opposed to the theoretical maximum value
in the flow. Last, we point out that Fig. 3 demonstrates the effect of varying d within the modified power-law model and of varying k for the Carreau model. It can be seen that the influence of d is effectively confined to regions of low shear whereas k affects the viscosity at all shear rates. The conclusion is that for order one dimensionless shear rate the modified power-law viscosity is essentially determined by the single parameter n whereas the Carreau viscosity is sensitive to the parameter k as well. We remark that, to the extent that these viscosity models provide an accurate representation of fluid rheology, ^ d and ^ k should be assumed fixed for a given fluid. Thus specifying either d or k in a numerical study risks being at odds with a laboratory realisation since ^ d and ^ k are then dependent on V and so will vary with the Reynolds number. However since d influences the viscosity only over a relatively limited range of shear rates we suggest that fixing d (rather than ^ d) in a numerical investigation should be of less concern than fixing k. Finally, we mention one potential disadvantage of adopting the modified power-law approach. Its use does require the introduction of the parameter j with non-physical units and complicates the definition of the Reynolds number since
Re ¼
V 2n an
j
:
ð7Þ
Nevertheless, as the widely-used power-law model has exactly the same drawbacks, we do not consider this to be a major disadvantage here.
128
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
3. Linear theory In order to validate our finite-element computations against results for an infinite aspect ratio channel we first revisit the wellknown system of ordinary differential equations describing the linear stability of a normal mode disturbance to axisymmetric flow driven by an azimuthal pressure gradient (see, for example, Drazin ðxÞ subject to a neutral perturbation and Reid [9]). For a basic flow v of the form
Þ þ eðu ~ ; v~ ; w; ~ p ~ÞE; ðu; v ; w; pÞ ¼ ð0; v ; 0; p ~ , etc. are functions of x only, the Navier–Stokes where E eiaz and u equations can be reduced to
a2 ðDDþ þ a2 Þ½l0 ðDDþ þ a2 ÞU 2Dðl0 Dþ UÞ 2Dþ ðl0 DUÞ 2l0 ~ u bv V; ¼ 2Dn2 r þ
Fig. 5. Neutral curves for the narrow gap Dean problem with d = 102 and n as shown.
r
Dþ ðlt D VÞ þ
lt r
ð8aÞ 2
D V a
l0 V ¼ ðDþ v ÞU;
ð8bÞ
~ ; v~ Þ ¼ ðU; ReVÞ. Here we have introduced the differential where ðu operator notation
D¼
d ; dx
Dþ ¼
d 1 d 1 þ ; D ¼ ; dx r dx r
together with the nondimensional parameter, the Dean number, defined as
Dn2 ¼ Re2 b1 : Within Eq. (8b) the quantities l0 and lt are, respectively, the viscosity associated with the basic flow and the tangent viscosity: n1
n3
l0 ¼ ð½D v 2 þ dÞ 2 ; lt ¼ ðn½D v 2 þ dÞð½D v 2 þ dÞ 2 : Note that if the basic flow viscosity were constant (as is the case for circular Couette flow in the narrow-gap limit) the viscosity terms could be taken outside the derivatives. For d = 0 (so that a true power-law model is recovered) and in the limit b ? 1 the first lt term in (8b) becomes the source of the additional factor of n found by Miller and Goddard [1]. Eqs. (8a) and (8b) together with boundary conditions
1 U ¼ DU ¼ V ¼ 0 on x ¼ ; 2
ð9Þ
constitute the linear stability problem for Dean vortices in a channel with a finite radius of curvature. The narrow gap linear stability problem is recovered by taking the limit b ? 1 such that D, D+ ? D and discarding all terms involving inverse powers of r in
(8a) and (8b) excepting the ratio b/r on the right hand side of (8a) which is set to one. The system (8a), (8b) and (9) is easily solved both for finite radii of curvature and in the narrow-gap limit. We do so using a shooting method based on a fourth order Runge–Kutta scheme and, for a fixed wavenumber a, obtain a critical Dean number, Dnc, for which the boundary conditions (9) can be satisfied. On varying a we obtain the function Dnc(a) describing the location of neutral stability. 3.1. Basic flow profile Before outlining the linear stability results we discuss briefly which, in the narrowthe associated non-Newtonian base flow v gap limit, is identical to the velocity profile of an equivalent plane Poiseuille flow. Fig. 4a shows the basic radial flow profile for several values of n. As the degree of shear thinning increases so the flow profile in the central portion of the channel flattens; moreover varying the value of d does not lead to a significant change in the solution. By way of illustration for all values of n in Fig. 4a the flow profile is shown for both d = 102 and 104. In each case the two are found to overlap to within the line width used. Fig. 4b shows the effective viscosity for the same values of n. As n decreases the central portion of the flow which is at very low shear, and so large effective viscosity, widens. Altering the value of d does, however, have an appreciable effect on the viscosity in this region since it determines its maximum value. 3.2. The narrow-gap case Fig. 5 shows the neutral stability curves obtained in the limit of narrow gap width, b ? 1, for d = 102 and selected values of n. For
for d = 102 and 104. (b) Effective viscosity l for d = 102. Fig. 4. Narrow gap baseflow for n = 0.2, 0.4, 0.6, 0.8 and 1. (a) Azimuthal velocity v
129
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
each value of n 6 1 the associated neutral curve has a minimum critical Dean number, Dnc(ac) at a critical wavenumber ac(n) in the range 3 < ac < 4. As n is decreased from unity this minimum Dnc(ac(n)) initially also decreases until n is approximately 0.5 at which point Dnc starts to increase again. This behaviour is in contrast to results for the Taylor problem for which both Dnc and ac are found to decrease monotonically with the shear-thinning index n (see Miller and Goddard [1]). However, we note that in the current formulation decreasing n does not lead to an everywhere reduced viscosity as we have nondimensionalised on the effective viscosity for c = 1. In fact decreasing n lowers the viscosity where c > 1 and increases it where c < 1. Since regions of both type are present in all flows considered here it is not enough to claim that increasing the intensity of shear thinning leads to a universally lower viscosity. Fig. 6 illustrates the effect of varying d on the stability characteristics of the flow. For values of n close to unity the effect of d is quite small. Once n 6 0.6 however the situation changes somewhat and the effect of varying d is much more significant, as can be seen in Fig. 6b where varying d over two orders of magnitude, which for n = 0.4 leads to a change in the peak viscosity by a factor of four, alters Dnc(ac) by well over 10%. This result is both surprising, given that d has a negligible impact on the basic flow, and inconvenient, in that it suggests that the stability properties of shear-thinning fluids are sensitively dependent on the precise constitutive relations that they satisfy and that a single governing index is not enough to capture this behaviour. It also raises the possibility that two fluids exhibiting apparently identical flow rates and rheology may have noticeably different stability characteristics due to their behaviour at shear rates too small to be measured without specialist equipment. Given the similarity of the basic flow it might be expected that the current findings will also apply to the stability characteristics
(a)
of plane Poiseuille flow. Nouar et al. [8] and Nouar and Frigaard [10] both examined this problem for a fluid obeying the Carreau model. However, these studies only discussed the effect of k upon the strength of the shear-thinning and made no attempt to separate the effect of varying k on the low-shear Newtonian plateau. Nouar et al. [8] did, however, highlight the significance of the Reynolds number definition in interpreting the effect on stability of shear thinning. On expressing their Reynolds number in the traditional manner, that is based upon the maximum attainable viscos^ 0 , they found that increased shear thinning (decreased n or ity l increased k) led to a less stable fluid. In contrast, choosing the tangent viscosity at the wall (the minimum viscosity present anywhere in the flow) as the reference viscosity implied that shear thinning actually stabilised the flow. It should be noted that this second definition leads to a Reynolds number that depends on k. This result indicates that care must be taken when drawing conclusions regarding the effect of shear thinning upon the stability of a flow and the associated physical consequences. It is apparent from Fig. 5 that the neutral curves for n = 0.2 and n = 0.4 contain inflection points not present for larger values of n. Fig. 6b shows that this effect becomes more pronounced as d decreases. In order to investigate this change in behaviour Fig. 7 shows the neutral eigenfunctions varying with a for the two values n = 0.6 and n = 0.2. In both cases there appears to be a transition with increasing wavenumber from an eigenfunction that fills the channel to one that is confined to x > 0. This is in line with the results of Hall [11] who demonstrated that for large wavenumbers the eigenfunction is located in an asymptotic region around a critical point xc lying in the outer half of the channel. A detailed discussion of the limit a ? 1 follows in Section 3.3. For n = 0.6 the eigenfunctions appear to vary relatively smoothly as a is increased. For n = 0.2 the flow is characterised by a much larger region of very high viscosity in the centre of
(b)
65 60 55
Dn
50 45 40 35 30 2
3
4
5
α
6
7
8
9
10
2
3
4
5
α
6
7
8
9
10
Fig. 6. Effect of varying d on the neutral curves for the narrow gap width Dean problem. Results are for d = 102 (solid line), d = 103 (dashed line) and d = 104 (dotted line). (a) n = 0.8. (b) n = 0.4.
Fig. 7. Variation with a of the neutrally stable eigenfunction for the narrow gap Dean problem. (a) n = 0.6, d = 104. (b) n = 0.2, d = 104.
130
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
the channel and the transition from an eigenfunction that spans this region to one that is contained in the outer portion of the channel appears to be much more abrupt. Decreasing n increases the size of this high viscosity region whilst decreasing either n or d increases the maximum viscosity there. Both actions cause the inflection point in the neutral curve to become more pronounced and it appears that it is the rapid transition of the eigenfunction across the high viscosity central portion of the channel that is largely responsible for the unusual shape of the neutral curves reported here.
l0 ¼ jV 1 jn1 ; l1 ¼ 2l0 ðn 1ÞV 1 1 V 2; 2 2 l2 ¼ 2l0 ðn 1Þðn 2ÞV 1 V 2 þ 3l0 ðn 1ÞV 1 1 V 3: We take expansions for the perturbation as follows: 1
1
V ¼ e2 ðV 0 ðgÞ þ e2 V 1 ðgÞ þ Þ;
U ¼ U 0 ðgÞ þ e2 U 1 ðgÞ þ ;
and for the governing parameter:
1 2Dn2 ¼ e4 d0 þ e2 d1 þ :
ð13Þ
Substituting these expansions into (11) we obtain, 3.3. Large wavenumber asymptotics
1
d0 ¼ l20 ðV 0 V 1 Þ ;
Hall [11] used a large wavenumber asymptotic solution to obtain an approximation to the right hand branch of the neutral curve for both the Dean and the Görtler problems. For the Dean problem with critical parameter equal to k0 at leading order he established that the WKBJ method was valid only in regions such that 1 k0 > ðv v 0 Þ ;
ð10Þ
where a prime indicates differentiation with respect to x. When v 0 Þ1 k0 ¼ ðv min the WKBJ roots coalesce at a second-order turning point located at a critical point xc and the method breaks down. Hall [11] noted that xc is the location at which the flow violates Ray 0 j is a maxileigh’s instability criterion [9] most strongly since jv v mum there. He demonstrated that in this case solutions to the disturbance equations can still be constructed in terms of parabolic-cylinder functions in a region of thickness e1/2 about the critical point
for e = a1 a small parameter. Furthermore these solutions represent the most dangerous modes available. We now proceed to apply the method of Hall [11] to the nonNewtonian Dean problem. For simplicity we assume a narrow gap approximation and a pure power-law fluid (d = 0), since the region in which we consider expansions will be removed from the low shear plateau. In this regime the linear stability equations (8a) and (8b) reduce to:
ð15Þ
at leading order. An application of (12) at first order then yields the consistency condition d1 = 0. At second order we obtain the consistency condition:
ðn þ 2ÞD2 V 0 þ g2 3ð3 2nÞ
ð11aÞ
V1
þ ð5 2nÞ
V2
dv 2 l ; dx
and so, noting that l ¼ jv 0 jn1 , to a critical location for the instabil 0 jv 0 j22n is a minimum: ity such that v v
V0
d2 ¼ d0 ðn þ 2Þ 3ð3 2nÞ
V3 V1
þ ð5 2nÞ
V2 V0
Dn ¼ a2
1 ðd0 þ a1 d2 Þ 2
ð16Þ
:
1=2 ð17Þ
;
where d0 and d2 are obtained from the basic flow profile according to (14) and (16). 3.4. The finite-gap case Results for the properties of Dean vortices in channels with a finite radii of curvature are limited, even for Newtonian fluids. Finlay and Nandakumar [12] calculated both ac and Dnc(ac) to very high accuracy and provided some detailed comparison with the results
1000
100
n= 1
1
l ¼ l0 þ e2 gl1 þ eg2 l2 þ
0.8 0.6 0.4
where
and
!#12
Fig. 8 shows neutral curves for d = 104 and various values of n together with the two term asymptotic prediction for the neutral Dean number
v ¼ V 0 þ e gV 1 þ eg2 V 2 þ
x¼xc
d2 V 0 ¼ 0: d0
and for non-negative integer values of m. The most unstable mode is therefore:
1 2
10
;
V0 þ
!#12 " 1 V3 V2 d2 ¼ 2d0 þ m ðn þ 2Þ 3ð3 2nÞ þ ð5 2nÞ ; 2 V1 V0
ð12Þ
Note that the variable viscosity has severed the link between the critical point and the point at which Rayleigh’s criterion is most strongly violated. In order to expand the governing Eqs. (11) about the critical point we expand the basic flow and its viscosity in terms of g
j 1 d v Vj ¼ j! dxj
!
ð11bÞ
An application of the WKBJ method to these equations leads to roots
d 2 ðv v 0 jv 0 j22n Þ ¼ l2 ðð3 2nÞv v 00 þ ðv 0 Þ Þ ¼ 0: dx
V3
Solutions of this equation can be obtained in terms of the paraboliccylinder function and are such that they decay as jgj ? 1 when
Dn
e2 ðD2 þ e2 Þ½lðD2 þ e2 ÞU 4DðlDUÞ ¼ 2Dn2 v V nDðlDVÞ e2 lV ¼ ðDv ÞU: 2
U 0 ¼ l
"
1
g ¼ e2 ðx xc Þ;
ðK 2 þ 1Þ ðnK 2 þ 1Þ ¼ k0 v
ð14Þ
1 0V 1 V 0;
1
α
10
Fig. 8. Logarithmic plot of the neutral curves for various values of n and d = 104 together with the two term asymptotic prediction for the limit a ? 1 given by (17).
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
131
4. Finite aspect ratio: A numerical study 4.1. Finite element method
Fig. 9. Neutral curves varying with curvature ratio, b, for n = 0.6 and d = 102.
We now examine whether the infinite aspect ratio linear stability results of Section 3 are applicable to a domain of finite extent in the z-direction, as depicted in Fig. 1. Of particular interest is whether the stability characteristics remain sensitive to d when the base flow is no longer a function of radius alone and when disturbances of finite amplitude are admitted. The nondimensionalised governing equations for axisymmetric flow were given in Section 2.1. Here we solve (1a)–(1c) and (2) using a finite-element method. We adopt the same approach as in Haines et al. [14] (herein referred to as HDB) and so only outline the main points of the procedure. The system of equations is solved on the rectangular domain
Fig. 10. Dnc against n for a = 3.93, b = 19.5 and d = 102 (solid line) and d = 104 (dashed line).
of previous studies. Our results for n = 1 agree with those of Finlay and Nandakumar [12] to at least four decimal places in both Dnc and ac. For n – 1 our linear stability analysis confirms that, as for the Newtonian case [13], the critical Dean number is significantly greater than that for the narrow-gap approximation, even for quite large values of b. Fig. 9 shows neutral curves for n = 0.6, d = 102 and a range of curvature ratios. We suggest that only for b > 100 can the neutral curve effectively be considered independent of the radius of curvature. Critical values of the Dean number against n for a = 2p/ 1.6 = 3.93 and b = 19.5 are shown in Fig. 10. We note that Fig. 10 does not indicate the onset of instability, since, for each value of n, it gives the critical value of the Dean number for a fixed a = 3.93 rather than the minimum taken over all wavenumbers Dnc(ac(n)). The value of a chosen is intended for later use in Section 4.2 where the results for an infinite aspect ratio presented here will be compared with computations in a channel of finite aspect ratio C = 1.6. For the values of n considered the wavenumber a = 3.93 is, however, always of a similar magnitude to the critical values which lie in the range 3 < ac(n) < 4. Fig. 10 indicates that d has a significant impact on the stability of the flow for values of n away from one. It also serves to illustrate the key findings of our linear stability study which we summarise as follows:
1 1 6x6 ; 2 2
1 6 Z 6 1;
ð18Þ
where we have defined a new coordinate Z = 2z/C. This has the advantage of introducing C alongside b as an explicit parameter in the equations so that continuation in either is possible. The computational domain is sketched in Fig. 11 and the weak form used to represent 1a,1b,1c and 2 can be found in Appendix A. The computational mesh is uniformly spaced in the coordinates x and Z and comprises (initially) Nx elements in the radial direction and NZ elements in the axial direction. A typical mesh resolution was of the order Nx = 20, NZ = 40. Experiments were conducted with various mesh dimensions and it was checked that the results discussed below were obtained on a sufficiently refined grid. Moreover, if any particular calculations required finer resolution, it was possible, as necessary, to adaptively refine the mesh as part of the solution procedure. The Z2 error estimator developed by Zienkiewicz and Zhu [15] with a flux based upon the shear stress was used to control the size of the mesh elements. The finite-element library oomph-lib Heil and Hazel [16] was used to assemble and solve the system of algebraic equations associated with a weak form of the axisymmetric Navier–Stokes equations. The Q2P1 finite element (see Gresho and Sani [17], p. 554) was used throughout. 4.2. Single wavelength validation In order to validate our finite-element computations against those for an infinite aspect ratio channel we first considered the case of a channel subject to symmetry conditions
w¼
@u @ v ¼ ¼ 0; @Z @Z
Decreasing n from one initially lowers the neutral curve, and in particular decreases Dnc(ac). This trend is reversed for n [ 0.5. For n [ 0.8 the value of d significantly affects Dnc. To the best of our knowledge these findings remain valid regardless of the curvature ratio.
Fig. 11. The axisymmetric computational geometry in terms of the coordinates given in (18).
132
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
at Z = ± 1. These boundary conditions allow the computation of a single vortex pair, the emergence of which with increasing Dean number can be compared with that predicted by the infinite aspect ratio linear stability analysis. On the side walls, x = ± 1/2, all velocities are set to zero with the exception of a small azimuthal velocity (everywhere less than 105 times the peak azimuthal velocity) on the inner cylinder boundary. This serves to slightly disconnect the bifurcation at Dnc by preferentially selecting the branch of vortices with outflow along the centreline of the domain. Fig. 12 shows results for C = 1.6, b = 19.5 and selected values of n 6 1. The radial velocity at the centre of the domain u(0, 0) is used as a suitable measure of the solutions. The vertical lines in Fig. 12 indicate the position at which a loss of stability (to the specified wavenumber) is predicted by linear theory. The agreement between linear theory and the two-dimensional numerical results is found to be very good. In particular Fig. 12 confirms that the Dean number for the onset of a single vortex pair does not decrease monotonically with n. Fig. 13 reproduces Fig. 12 with additional data for d = 104. The agreement with linear theory is again found to be good, with the results in Fig. 13 confirming the increasing influence of d with decreasing n. Fig. 13 also exhibits some striking features not captured by the linear theory. Firstly, although the value of d significantly affects
Fig. 12. Solution branches for a single vortex pair of wavelength 1.6 contained in a finite domain subject to symmetry conditions at Z = ± 1. Results are for C = 1.6, b = 19.5 and d = 102. Vertical lines and accompanying values indicate the position at which a loss of stability is predicted by linear theory. The measure of the flow used is the centreline radial velocity at the point x = Z = 0.
Fig. 13. Solution branches for a single vortex pair of wavelength 1.6 contained in a finite domain subject to symmetry conditions at Z = ± 1. Results are for C = 1.6, b = 19.5 and d = 104. Results for d = 102, previously given in Fig. 12, are shown in grey.
the location of Dnc, it has less of an impact on the vortex modes once they become nonlinear. For all three values of n < 1 shown in Fig. 13 the solution branches for the two values of d move closer together as the strength of the vortices (indicated by ju(0, 0)j) increases. This can be attributed to a reduction in the amount of fluid at very low shear and to the corresponding decrease in importance of the details of the zero-shear plateau. The second striking feature of Fig. 13 is the presence of hysteresis for d = 104. This hysteresis becomes more pronounced as n is reduced. For n = 0.4 plausibly stable nonlinear vortex modes exist for Dean numbers in the range 35.2 < Dn < Dnc = 35.83. 4.3. Comparison with Shanthini and Nandakumar Shanthini and Nandakumar [2] explored the effect of the power-law index n on the bifurcation structure in a square duct. They performed symmetric computations that identified the stable two and four-cell parts of the main solution branch and noted that the non-Newtonian characteristics of the fluid did not appear to change the qualitative behaviour of the flow. The turning point at which the two-vortex solution jumps to a four-vortex solution with increasing Dn (Dn24 in Shanthini’s notation) and the corresponding point at which the four-vortex solution jumps to a two-vortex solution with decreasing Dn (Dn42) are, for n = 1, the bifurcation points L1 and L2 described by Winters [3] (see his Fig. 3). A recreation of Fig. 6 of Shanthini and Nandakumar [2] showing the location of the two bifurcation points in the Dn–n plane is given in Fig. 14. Our results, obtained by using bifurcation detection to track the turning points Dn24 and Dn42 over a range of n, are shown as solid and dashed lines respectively. The original data points of Shanthini and Nandakumar [2], computed to an accuracy of Dnc ± 1 using the bisection method, have been estimated graphically. The agreement is generally good, although there appears to be a departure for Dn24 when n = 0.6. At such small aspect ratios the flow experiences non-zero shear everywhere apart from exactly in the corners of the domain. Shanthini and Nandakumar [2] dealt with this corner singularity by imposing a maximum viscosity of 106 there, far greater than was found anywhere else in the flow. Although they do not appear to have performed any sensitivity analysis on their viscosity ceiling we have observed that for a square duct and d 6 104 the correction to a purely power-law viscosity is felt only in the corners and that any further reduction in d does not influence the bifurcation loci significantly. Thus the parameter d is found to be of
Fig. 14. Loci of fold bifurcations in the Dn–n plane for b = 100 and C = 1. Solid line: Turning point at which two vortex solution jumps to four vortex solution with increasing Dn, Dn24 in Shanthini’s notation. Dashed line: Turning point at which four vortex solution jumps to two vortex solution with decreasing Dn, Dn42 in Shanthini’s notation. Shanthini’s data points are also shown, they have been reconstructed from Fig. 6 of the original paper.
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
133
limited importance for such small aspect ratios. Consequently the power-law model can be, and has been, used without the care that we shall see is required for larger aspect ratios, for which regions of very low shear persist in the interior of the domain.
5. Large aspect ratio results We now consider the axisymmetric solution structure for curved channels of large aspect ratio, which we take to mean C 10. The corresponding structure for a Newtonian fluid was previously described in HDB. For C = 10 we solve the governing equations for flow in the upper half channel Z > 0 subject to the symmetry boundary conditions w = uZ = vZ = 0 on Z = 0. The velocities on all other boundaries are set to zero. We first compute the unique low Dean number solution and then use arc-length continuation to track the associated (primary) solution branch to higher Dean number. Fig. 15 depicts this solution branch close to the onset of instability for a modified power-law fluid with shear-thinning index n = 0.6 and d = 102. Instability is found to occur at a lower Dean number than for the Newtonian case, as suggested by the linear stability calculations in Section 3. Fig. 15 also shows streamline plots of the secondary velocity at selected points. Note that in all cases the flow is symmetric about Z = 0 and so only the half of the domain Z > 0 is shown. As for a Newtonian fluid the flow at low Dean number is characterised by a large end vortex but an absence of vortex motion on the channel centreline (u(0, 0) 0). As the Dean number is increased through the critical value for centrifugal instability the flow on the primary branch changes to one with vortices throughout the channel. Fig. 15 also illustrates two disconnected solution branches. We refer to the solutions on these branches as secondary modes. A streamline plot shows the crosssectional flow for one of these modes and it can be seen that the number of vortices in the half channel differs by one from the primary mode. The (secondary) branch of solutions characterised by outflow along the centreline is formed of two distinct parts. Whilst this was not the case for a Newtonian fluid at the same aspect ratio, similar behaviour was observed in the Newtonian case at slightly smaller aspect ratios (see HDB). We shall return to this observation in Section 5.1 below. The effect of decreasing d upon the bifurcation structure depicted in Fig. 15 is shown in Fig. 16. Reducing d by two orders of magnitude (corresponding to an increase in the fluid’s zero-shear
Fig. 16. Symmetric solution branches for n = 0.6 and for d = 102 (solid line), 103 (dashed line) and 104 (dotted line). For d = 103 and 104 we have made no attempt to continue the solution branches beyond Dn 29.4.
viscosity of roughly 2.5) leads to a change of around 4% in the critical Dean number. This is qualitatively in line with both the linear stability results of Fig. 10 and the single wavelength results of Section 4.2. We do not expect quantitative agreement for a number of reasons discussed in HDB. These include the difficulty of establishing a wavenumber for the onset of the finite aspect ratio instability and a discrepancy in the magnitude of the nondimensionalised centreline velocity between the finite and infinite aspect ratio results. Both of these phenomena may be attributed to the presence of end effects. Although an aspect ratio of 10 may be regarded as large it is clear from Fig. 15 that the end vortex region still takes up a significant portion of the channel both with and without the presence of the instability. As in Section 4.2 the results in Fig. 16 also enable the effect of d on nonlinear perturbation growth to be examined. As was noted for a single wavelength, reducing d introduces a hysteresis region to the primary solution branch. The size of this hysteresis region is observed to increase with decreasing d and also with decreasing n (although this effect is not shown here). Again, just as in Section 4.2, we find that the impact of d is at its most pronounced close to the onset of instability. As the strength of the vortices grow throughout the channel (indicated by increasing ju(0, 0)j) d has a decreasing impact on the solution branches due to a reduction in size of the low shear regions where its effect is felt.
Fig. 15. Symmetric solution branches for a modified power-law fluid with n = 0.6 and d = 102 and for C = 10, b = 19.5. The labels attached to the limit points correspond to those in Fig. 17. The inset streamline plots show the secondary flow found on selected solution branches.
134
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
Fig. 17. The location of limit points in the Dn C plane for (a) n = 0.6 and d = 102 and (b) the Newtonian fluid as reported previously in HDB. Note that the two figures cover different coordinate ranges. The labels I and II relate the sets of limit points shown in (a) to those given (for C = 10) in Fig. 15.
5.1. Bifurcation structure with varying aspect ratio Neither the solution branches or secondary streamlines depicted in Fig. 15 bear much resemblance to those for a Newtonian fluid at the same aspect ratio (see, for example, Fig. 3 of HDB). However, the rapid change of these features with aspect ratio observed by the previous study indicates that a description of the flow at a single aspect ratio does not tell the full story. Hence now we detail results for the loci of limit points in the Dn C plane. Corresponding loci for the Newtonian problem were given in HDB where it was found that the primary mode of instability (characterised by the number of vortices accommodated in the channel cross-section) was selected via a transcritical bifurcation between two different modes. This process was first described by Benjamin [18] in the context of the Taylor problem. In order to examine the non-Newtonian bifurcation structure the two turning point bifurcations shown in Fig. 15 (for n = 0.6 and d = 102) have been continued over a range of aspect ratios 6 < C < 13. The results are shown in Fig. 17a. A comparison with Fig. 6 of HDB (reproduced in part here as Fig. 17b) demonstrates that qualitatively similar features are found at an aspect ratio somewhere between one and two greater for n = 0.6 and d = 102 than in the Newtonian case. For example the bifurcation structure for C = 10 given in Fig. 15 is found to be similar to that found for a Newtonian fluid with C 8 (see the middle schematic of Fig. 7a of HDB). Although several other turning point bifurcations were found for the Newtonian problem we have made no attempt to compute them here. As Shanthini and Nandakumar [2] observed for C = 1, shear thinning does not appear to change the qualitative structure of the bifurcation set. 6. Discussion We introduced a new viscosity function (3) with similarities to both the power-law and Carreau models. In addition to the shearthinning index, n, this viscosity function contains a second parameter d that controls the location of the low-shear Newtonian plateau. This second parameter has the advantage over the Carreau model that it does not significantly affect the viscosity at dimensionless strain rates of order one. Nevertheless, and despite having a negligible impact on the basic flow profile d is found to have a significant impact upon the onset of instability. This linear stability finding is confirmed by fully nonlinear results for axisymmetric flow in a finite aspect ratio duct, computations that additionally demonstrate the existence of hysteresis on the primary branch of vortices. Decreasing either n or d increases the size of the
hysteresis region. As the nonlinear vortices increase in strength so the proportion of the domain at very low shear decreases and consequently so does the influence of d. Although the structure in Fig. 15 appears to be very different from that for a Newtonian fluid (Fig. 3 of HDB) a bifurcation analysis shows that this can be attributed to a shift in the aspect ratio at which qualitatively similar features are present rather than to a fundamental change in the character of the vortices. Given the sensitive dependence upon aspect ratio of the vortex structures present in both the Dean and Taylor–Couette problems, bifurcation tracking provides a powerful tool with which to unravel the underlying behaviour. Finally we speculate as to whether plane Poiseuille flow exhibits a similar sensitivity to the low-shear Newtonian plateau. Although the nature of the two instabilities is fundamentally different the basic flow profile in the two problems is very similar being characterised by a velocity profile that is a function of a single transverse coordinate and that has a point of vanishing shear somewhere in the interior of the flow. Whilst the instability of plane Poiseuille flow of a shear-thinning fluid has been previously studied [8,10] the use of the Carreau model meant it was not possible to isolate the effect of k upon the Newtonian plateau from its effect on the viscosity at larger shear rates. Acknowledgements JPD and APB gratefully acknowledge the support of the Australian Research Council. The comments of the anonymous referees are gratefully acknowledged. Appendix A. Weak form In the finite-element computations Eqs. (1a)–(1c) and (2) were represented by the following weak forms:
Rui ¼
Z Z
@u @u w @u v 2 þu þ Wi r dx dZ @t @x r @Z r X Z Z p u 2l 2 Wi r dx dZ r X r Z Z Z Z @u @ Wi 1 @u @w 1 p 2l l þ r dx dZ þ @x @x r @Z @x r X X @ Wi
r dx dZ; @Z St
where we define Rui ð1 6 i 6 NÞ to be the residual associated with evaluating the r-momentum equation with respect to the test function Wi, r = C/2 and X is our computational domain (18). Similarly:
P.E. Haines et al. / Journal of Non-Newtonian Fluid Mechanics 198 (2013) 125–135
Z Z @v @ v w @ v uv Rvi ¼ St þu þ þ Wi r dx dZ @t @x r @Z r X Z Z l @v v 1 @p Wi r dx dZ r @h r @x r X Z Z Z Z @ v v @ Wi 1 @v 1 þ l l r dx dZ þ @x r @x r @Z r X X @ Wi
r dx dZ; @Z and
Rw i ¼
Z Z @w @w w @w St þu þ Wi r dx dZ @t @x r @Z X Z Z Z Z @w 1 @u @ Wi l @w 1 l þ p2 þ r dx dZ @x r @Z @x r @Z r X X @ Wi
r dx dZ; @Z
The weak form for conservation of mass is constructed using the (linear) pressure basis functions Uj (1 6 j 6 M) as test functions:
Rpj ¼
Z Z @u u 1 @w þ þ Uj r dx dZ: r r @Z X @r
References [1] C. Miller, J. Goddard, Appendix to polymer fluid mechanics, Adv. Appl. Mech. 19 (1979) 143–219. [2] W. Shanthini, K. Nandakumar, Bifurcation phenomena of generalized Newtonian fluids in curved rectangular ducts, J. Non-Newton. Fluid Mech. 22 (1) (1986) 35–60.
135
[3] K. Winters, Bifurcation study of laminar flow in a curved tube of rectangular cross-section, J. Fluid Mech. 180 (1987) 343–369. [4] H. Fellouah, C. Castelain, A. Ould-El-Moctar, H. Peerhossaini, A numerical study of Dean instability in non-Newtonian fluids, J. Fluids Eng. 128 (2006) 34–41. [5] F. Delplace, J.C. Leuliet, Generalized Reynolds number for the flow of power law fluids in cylindrical ducts of arbitrary cross-section, Chem. Eng. J. Biochem. Eng. J. 56 (1995) 33–37. [6] H. Fellouah, C. Castelain, A. Ould-El-Moctar, H. Peerhossaini, The Dean instability in power-law and Bingham fluids in a curved rectangular duct, J. Non-Newton. Fluid Mech. 165 (3–4) (2010) 163–173. [7] R. Bird, R. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, Fluid Mechanics, vol. 1, John Wiley and Sons Inc., New York, NY, 1987. [8] C. Nouar, A. Bottaro, J. Brancher, Delaying transition to turbulence in channel flow: revisiting the stability of shear-thinning fluids, J. Fluid Mech. 592 (2007) 177–194. [9] P.G. Drazin, W.H. Reid, Hydrodynamic Stability, Cambridge University Press, 1981. [10] C. Nouar, I. Frigaard, Stability of plane Couette–Poiseuille flow of shearthinning fluid, Phys. Fluids 21 (2009) 064104. [11] P. Hall, Taylor–Görtler vortices in fully developed or boundary-layer flows: linear theory, J. Fluid Mech. 124 (1) (1982) 475–494. [12] W. Finlay, K. Nandakumar, Onset of two-dimensional cellular flow in finite curved channels of large aspect ratio, Phys. Fluids A: Fluid Dynam. 2 (1990) 1163–1174. [13] E. Sparrow, On the onset of flow instability in a curved channel of arbitrary height, Z. Angew. Math. Phys. (ZAMP) 15 (6) (1964) 638–642. [14] P.E. Haines, J.P. Denier, A.P. Bassom, The Dean instability for finite aspect ratio ducts, J. Fluid Mech. 716 (2013). R8-1-13. [15] O. Zienkiewicz, J. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique, Int. J. Numer. Methods Eng. 33 (7) (1992) 1331–1364. [16] M. Heil, A. Hazel, oomph-lib-An object-oriented multi-physics finite-element library, in: M. Schafer, H.-J. Bungartz (Eds.), Fluid-Structure Interaction, Lecture Notes on Computational Science and Engineering, 2006, pp. 19–49. [17] P. Gresho, R. Sani, Incompressible Flow and the Finite Element Method, Advection-Diffusion and Isothermal Laminar Flow, vol. 1, John Wiley & Sons, 1998. [18] T. Benjamin, Bifurcation phenomena in steady flows of a viscous fluid. I. Theory, Proc. Roy. Soc. Lond. A. Math. Phys. Sci. 359 (1696) (1978) 1–26.