Fully three-dimensional, time-dependent numerical simulations of Newtonian and viscoelastic swirling flows in a confined cylinder

J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

Fully three-dimensional, time-dependent numerical simulations of Newtonian and viscoelastic swirling flows in a confined cylinder Part I. Method and steady flows 夽

S.C. Xue, N. Phan-Thien ∗ , R.I. Tanner Department of Mechanical and Mechatronic Engineering, The University of Sydney, Sydney, NSW 2006, Australia Received 18 June 1999; received in revised form 14 July 1999

Abstract Aiming at simulating the viscoelastic flow problems involving three-dimensional, time-dependent complex flow, a fully three-dimensional (3D) numerical method using an implicit finite volume formulation is developed in the cylindrical coordinate system. We focus our attention on the well known swirling flows in a confined cylinder with a rotating bottom lid, and through careful comparisons with the available experimental results, the suitability of different constitutive models used to describe viscoelastic fluids is evaluated. Three difficulties arising from the coordinate transformation in a structured mesh system are removed by some simple, but effective measures. With this method, the velocity fields measured in experiments are accurately predicted, and the development of vortex breakdown in terms of Reynolds number (Re) and the height to radius ratio (H/R), as well as the time, observed in experimental visualisations are exactly reproduced for Newtonian flows. The influence of elasticity on the vertex breakdown processes is investigated for the low viscosity dilute polyacrylamide (PAA) solutions described by the upper converted Maxwell (UCM) model. It is confirmed numerically that the existence domain of vortex breakdown occurs at significantly greater Re than that for Newtonian solvent of similar shear viscosity. It seems that the centrifugal force leading to vortex breakdown in Newtonian fluids is balanced by the normal stresses in elastic fluids. As a result, vortex breakdown would be significantly delayed in terms of Re or even suppressed with increasing elasticities. Numerical simulations of the steady swirling flow involving viscoelastic fluids with a high constant viscosity (low Re) shows that the flow field could partially or entirely reversed depending on the level of fluid elasticity. Our numerical results demonstrate that at a low value of Re, for the fluids with a medium high elasticity characterized by the elastic number El = O(1), an elastic ‘ring’ vortex is developing and growing with increasing elasticities at the edge of the rotating lid, and trying to push the Newtonian vortex to the radial centre of the rotating lid, and finally it occupies the whole flow domain at a certain level of elasticity. It is also found that a small weak edge ring vortex appears at the edge of the rigid cover and counter-rotates with the main vortex when the aspect ratio H/R ≥ 1, and its size tends to increase with increasing H/R. 夽

Dedicated to Professor David V. Boger on the occasion of his 60th birthday. *Corresponding author. Tel.: +61-2-9351-7127; fax: +61-2-9351-7060 E-mail address: [email protected] (N. Phan-Thien) 0377-0257/99/$ – see front matter ©1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 0 2 5 7 ( 9 9 ) 0 0 0 7 3 - 7

338

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

The double-cell flow structure for strong shear-thinning viscoelastic fluids observed in visualisation has been successfully predicted numerically by using the simplified Phen-Thien–Tanner (SPTT) model. With steady-state calculations, at a given inertia level, when the flow field is dominated by elasticity, we found that there is a small counter-rotating central ring vortex appears around the centre of the rotating lid, which was usually reported as a sign of instability of the flow field in experimental work. In the second part of this work, we will focus on the numerical inspection of the highly unsteady spiral vortex flows of viscoelastic fluids observed in experiments by using time-dependent calculations, thus, checking the stability criterion of McKinley et al. [G.H. Mckinley, P. Pakdeland, A. Öztekin, J. Non-Newtonian Fluid Mech. 67 (1996) 19–47], and evaluating the suitability of different constitutive equations. ©1999 Elsevier Science B.V. All rights reserved. Keywords: Polyacrylamide solutions; Upper converted Maxwell model; Newtonian fluids

1. Introduction Viscoelastic flow simulations represent a major challenging task in Computational Fluid Dynamics (CFD), especially in three-dimensional (3D) and transient flow problems. Obviously the efficiency in terms of memory requirement and CPU time consumption, as well as the reliability of the method are critical issues. Another issue related to the performance of the method is the so-called high Weissenberg number (We) problem [2]; that is, the calculations fail to converge at a moderate level of fluid elasticity. This is especially true when the flow involves some singularities, either geometrically or in the stress [3]. More importantly, the choice of the constitutive equations must be considered with care, if the simulation is to reproduce the viscoelastic behaviour observed in experimental work. The traditional way to evaluate a constitutive equation is to check its predictions in simple flows, for example, shear or elongational flow. However, in most flows of practical interest, both shear and elongational deformations are involved. Obviously, the information provided by fitting simple flows may not be enough to ensure reliable predictions in complex flows. Therefore, complex flows should be used to test the predictability of constitutive equations. To this end, it is desirable to set up a complex (non-viscometric) flow problem in which the boundary conditions can be defined unambiguously, and for which experimental measurements and visualisations can be easily carried out so that direct comparison between numerical predictions and experimental observations can be made. Confined swirling flow in an enclosed cylindrical vessel with a rotating lid, as shown in Fig. 1a, provides such a good test case [6]. This flow system consists of a fixed enclosed circular cylinder of radius R and height H with a lid of radius R fitted across the base. The cylinder is completely filled with liquid and the lid is impulsively set to rotate with a constant angular velocity . A non-uniform centrifugal force, as well as the normal stresses if the liquid is viscoelastic along the base is produced along with a tangential primary flow due to the rotation of the lid, thus causing a weak superimposed toroidal secondary flow whose direction of rotation depends upon the properties of fluid. After a transient phase from an initial rest state, depending on the geometric and physical parameters involved (the ratio of height and radius H/R, the Reynolds number Re, and rheological properties of fluid if it is viscoelastic), a steady-state and possibly time-dependent periodic flow situation may be reached. If the fluid is Newtonian, then it will be driven outwards along the lid by secondary flow, expelled up the side walls and drawn down the central axis to maintain continuity in a spiral motion. Under certain conditions, an undulating meridional flow along with one or several recirculatory regions (separation

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

339

Fig. 1. (a) Schematic diagram of an enclosed cylinder with a rotating lid in the cylindrical coordinate system (r,θ,z) and the mesh structure, and (b) a typical vortex breakdown pattern for a Newtonian fluid in H/R = 2.0 at Re = 1854.

bubbles) along the rotation axis appear. This flow pattern is interpreted as a ‘vortex breakdown’ which is typically shown in Fig. 1b in a cylinder with H/R = 2 at Re = 1854. The phenomenon of vortex breakdown in the confined swirling flow of Newtonian fluids has been widely observed in experimental visualisations (e.g., [7,14]) and investigated numerically or theoretically (e.g., [8–13]). Escudier [7] extended the experimental observations by Vogel [4] and Ronnenberg [5] and reported that, depending on the two dimensionless variables, Re and R/H, multiple breakdown bubbles may co-exist in the enclosed cylindrical geometry. They found that the recirculation bubbles are axisymmetric and steady over a large range of governing parameters. An existence domain for the occurrence of breakdown was depicted in terms of Re (1000–3500) versus R/H (1.2–3.7). By solving the Navier–Stokes equations using a finite-difference method with the stream function-vorticity formulation, as well as the assumption of axisymmetry, the existence domain has been faithfully reproduced numerically (see, [9–13]). However, in a recent experimental work, Spohn et al. [14] claimed that the breakdown bubbles are ‘open’, and highly axisymmetric on the upstream side of the bubble and asymmetric on their rear side, and they argued that their observations were real and was not caused by the visualisation technique. If that is the case, then numerical calculations of the confined swirling flow problem have to be carried out using a 3D numerical method without the assumption of axisymmetry. Much more complex flow patterns may be induced for non-Newtonian fluids as a result of their elasticity. The most striking difference is that the toroidal secondary flow may reverse, and become unstable even at low rotation rates depending on the fluid elastic properties. The first experimental investigation of the flow patterns of viscoelastic fluids in the enclosed cylinder system was due to Hill et al. [6], who observed the reversed toroidal secondary flow at low rotation rates for shear thinning elastic fluids (polyacrylamide in a solvent of glycerol amid water). For some combinations of Re and We, the secondary flow is a

340

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

complicated compromise between the two types of flows [15], that is, there are two counter-rotating vortices, each struggling to be dominant. A highly unsteady flow pattern was observed at relatively high rotation rates in which the flow changes direction in a confused manner, and a small vortex located at the centre of the rotating lid was spun off into the main body of the fluid every 20 s and subsequently reformed at where it was. Using an aqueous shear thinning elastic fluid, Day et al. [16] also observed flow reversal at low rotation rates, and as increasing rotation rates, an instability with the formation of a spiral vortex flowing away from the rotating lid was observed. Using strongly shear thinning, but slightly elastic fluids, Escudier and Cullen [17] investigated experimentally the influence of shear thinning on the viscoelastic swirling flow structures. At well below a critical value of Re (for any given H/R) at which vortex breakdown would be expected to occur, they observed the double-cell structure in which the lower cell close to the rotating lid is stronger and Newtonian-like, and the upper cell close to the enclosed top is weaker and reversed. Obviously, the occurrence of such a structure is the result of compound effects of inertia, shear thinning as well as elasticity. They also noticed that an upflowing jet of fluid close to the axis of symmetry is in evidence in many instances, which was thought as an indication of instability. As we know, the factors influencing the flow behaviour include elasticity and shear-rate-dependent viscosity. To investigate the exclusive effects caused solely by elasticity, thus eliminating the possible effects of shear thinning on the flow behaviour, a series of remarkable experimental work using Boger fluids [18] have been done by Stokes et al. ([19–21]). Using a polyacrylamide (PAA) Boger fluid with a high viscosity (to minimise the effects of fluid inertia), they [19] confirmed that the observed phenomena is caused purely by the elastic nature of the fluid. Using highly dilute flexible PAA solutions with a constant low viscosity such that inertia effects are not negligible and tend to dominate the flow, they [20] found that the flow is ‘Newtonian-hike’, but the occurrence of vortex breakdown observed in a Newtonian solution of a similar viscosity has been postponed, and even suppressed with increasing polymer concentrations. Recently, using a range of different fluids, they [21] observed that the swirling flow of PAA Boger fluids became unstable under the conditions ranging from being inertia-dominated to elasticity-dominated, and the values of Re required for the occurrence of the instabilities are correlated extremely well with the elasticity number El defined as the ratio between We and Re. They suggested that the instability is likely to be related to the elastic instabilities observed in other rotating flow, such as Taylor–Couette flow. Again, due to the 3D and unsteady nature of the spiral flow, a fully 3D simulation is preferred for the viscoelastic cases. However, numerical probing into the swirling flow problem involving viscoelastic fluids has been confined to axisymmetric flows and the use of relatively simple explicit models. They include Krammer and Jonson’s perturbation analysis [22] with the WJFLMB model, Nirschl and Stewart’s orthogonal collocation fully [23] and Chiao and Chang’s global spectral method [24] using the CEF model, and Wunsch and Bohme’s mixed finite element method using the Wagner model [25]. Recently, more realistic non-linear Giesekus model was first used to describe viscoelastic fluids in simulation of the confined swirling flow by Moroi et al. [26]. With simple models, only qualitative predictions are expected. Krammer and Jonson’s perturbation analysis [22] with the WJFLMB model indicated qualitatively that the primary normal stress function has an important influence on the steady flow field, while the second normal stress function is less significant. With increasing fluid elasticity, the flow field may partially reversed and finally entirely reversed compared to the corresponding Newtonian field. Using an orthogonal collocation method and the CEF model to fit the viscoelastic fluids used in Hill’s experiment, Nirschl and Stewart [23] was able to predict qualitatively some of the steady-state

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

341

observations. But, they noticed the possibility of a Hopf bifurcation in this system which lead to a time-periodic solution. Using a global spectral method, Chiao and Chang [24] extended Nirschl and Stewart’s analysis. For steady-state flows, they qualitatively mapped out Hill’s measurements of flow fields and stability boundaries for the different regions including counter-rotating vortices, a single Newtonian vortex and a single viscoelastic vortex in the rotation rate and elasticity concentration parameter space. Their nonlinear analysis demonstrated the existence of a region of chaotic flow, but they were not able to locate the time-periodic solution in the unstable region because of the limitation of the numerical method used. To our knowledge, no numerical studies on unsteady complex flows of viscoelastic fluids in this geometry have been yet carried out. This flow is simple in geometry and rich in its information about the flow fields, and with a well of experimental data already available, is most suitable to evaluate the suitability of a particular constitutive equation. Obviously, to capture such complex flows: 1. fully 3D calculations are required without the assumption of axisymmetry which precludes nonaxisymmetric instability; 2. the flow problem has to be taken as an initial value problem. That is, the fully Navier–Stokes equations and non-linear constitutive equations (for the viscoelastic cases) must be integrated in time to predict the developing processes of the flow fields from the start-up state to the stable region, and to locate the time-periodic solution in the unstable region. In recent years, increased research on numerical analysis of viscoelastic flows and the improved computer performance permit more efficient solutions of the complete 3D Navier-Stokes and constitutive equations. Of the many numerical methods in use, finite volume method (FVM) has made 3D viscoelastic flow problems more tractable with its time and space saving features. However even though this method has been extended to a large variety of two-dimensional (2D) viscoelastic flow problems and some 3D problems (see, Xue et al. [28,29]). it has hot been used for 3D viscoelastic problems involving cylindrical geometry. The aim of this paper, which will be presented in two parts, is to extend the efficient fully 3D implicit FVM for viscoelastic flow simulations developed previously to the problems involving cylindrical geometry, and carry out 3D, time-dependent calculations on the confined swirling flows of Newtonian and viscoelastic fluids without the assumption of axisymmetry. The typical effects of elasticity and share-rate-dependent viscosity are uncoupled in numerical simulations. By comparing the predicted results with the corresponding experimental measurements and visualisations, the efficiency and accuracy of the method, and the suitability of a particular constitutive model used to describe viscoelastic fluids will be fully evaluated. The first part of the paper presents the numerical method and its performance in terms of efficiency and accuracy by simulating the following flows and comparing with the available experimental results whenever possible: 1. the transient and steady-state vortex breakdown process for Newtonian fluid; 2. the steady flow of constant viscosity elastic fluid described by the UCM model [27] and 3. the steady flow of strong shear-thinning, low elastic fluid described by the simplified Phan-Thien–Tanner (PTT) model [30]. The unsteady viscoelastic flows will be presented in Part II. Note that the decision to use the UCM model is motivated by the fact that the UCM model is difficult to handle numerically — any numerical method that can handle the UCM model effectively should be able to do well with other more complex constitutive equations — and that the flow is mainly shear.

342

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

2. Problem formulation We wish to simulate the incompressible, isothermal and time-dependent complex 3D flow of viscoelastic fluid in a fixed enclosed circular cylinder of radius R and height H with a lid of radius R fit across the base. The cylinder is completely filled with fluid and the lid is impulsively set to rotate with a constant angular velocity . The geometry of the flow problem is sketched in Fig. 1a. 2.1. Governing equations The equations governing the flow can be written as follows: • Continuity equation (conservation of mass): ∇ ·u = 0

(1)

• Momentum equation (conservation of momentum in absence of body forces): ∂ u) + ∇ · ρuu uu = ∇ · (−pII + σ ) (ρu (2) ∂t where t is the time, ρ the mass density, and u the velocity vector. In this expression, the total stress is decomposed into two parts: the hydrostatic pressure with I being the unit tensor, which produces no stress power in an isochoric motion, and another part is called the extra stress tensor denoted by ␴ which is related to kinematic quantities by an appropriate constitutive equation. • Constitutive equation In this paper, we are mainly interested in differential models in which the total extra stress is arbitrarily split into two components: σ = 2ηN D + τ

(3)

For viscoelastic fluids, the stress τ is assumed to satisfy the PTT model: ∇

gττ + λτ = 2ηm0D ,

(4)

with ∇

τ =

∂ττ + u · ∇τ − L · τ − τ · LT , ∂t

(5)

and g =1+

λε tr τ ; ηm0

D, L = L − ξD

(6)

uT is the velocity gradient tensor, T denotes the transpose operation, λ the relaxation time, where L = ∇u ηm0 is the zero shear rate molecular-contributed viscosity. The zero-shear rate viscosity of the fluid is ηm = ηN + ηm0 . By defining β = ηm0 /η0 as the retardation ratio, we have ηN = (1 − β)η0 and ηm0 = βη0 . The constants and ξ , the model parameters. For zero ξ , the model is reduced to the simplified PTT (SPTT) model [31]. 1n the case when ε = 0 and β = 1, the UCM model is recovered. For viscoelastic calculations, from the numerical stability point of view, it is desirable to make the elliptic operator in the–momentum equation as large as possible [32]. For this purpose, various

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

343

schemes designed to treat the coupling of the kinematics and extra stresses have been developed, such as the elastic-viscous split stress (EVSS) formulation of Rajagopalan et al. [33], and the adaptive viscoelastic stress splitting (AVSS) scheme of Sun et al. [34]. In this work, the EVSS formulation is incorporated into our 3D FVM method [28]). With the EVSS formulation, the polymer-contribution stress ␶ was expressed in terms of the elastic stress tensor τije , that is, τ = τ e + 2Eβη0D where E=

(

1

with EVSS

0

without EVSS

(7)

(8)

By substituting this into the constitutive equation [Eq. (4)], we obtain the equation for ␶ e : ∂ uτ e ) = 2βη0 (1 − Eg)D D − gττ e + λ[L · (ττ e + Eττ N ) + (ττ e + Eττ N ) · LT ] (λττ e ) + ∇ · (λuτ ∂t   ∂ N N uτ ) , −E (λττ ) + ∇ · (λuτ ∂t

(9)

where τ N = 2βη0D ,

(10)

and the momentum equation becomes ∂ u) + ∇ · (ρuu uu − η0 ∇u u) = −∇pII + ∇ · [ττ e − 2(1 − E)βη0D ]. (ρu ∂t

(11)

Thus, an elliptic operator is kept in the momentum equation with the extra term now being ∇·␶e , and the dependent variables of the system are now u , p and ␶ e (in the rest of the paper, the superscript will be dropped for sake of simple expression). Adding an artificial diffusion term on both sides of the constitutive equation is another way to stabilise viscoelastic calculations. However, in view of its under-relaxation function (see, [28]), it is not suitable for time-dependent calculations. 2.2. Dimensional analysis and boundary conditions Experimentally, it has been shown that the features of the confined swirling flow are well characterised by several nondimensional numbers. To nondimensionalize the physical variables, it is convenient to introduce R and R as the characteristic length and velocity, thus, the time and stress being scaled by 1/ and η0  The momentum and constitutive equations in dimensionless form can be recast as   D] u −∇pII + ∇ · [τ − 2(1 − E)βD 1 ∂u u = ∇u (12) + ∇ · uu − ∂t Re Re and

344

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

2β(1 − Eg) − gττ ∂ττ uτ ) = + ∇ · (uτ + [L · (ττ + Eττ N ) + (ττ + Eττ N ) · LT ] ∂t We # " N ∂ττ N uτ ) , + ∇ · (uτ −E ∂t

(13)

where Re =

ρR 2  , η0

We = λ

(14)

are Reynolds number and Weissenberg number, respectively. Obviously, for Newtonian fluid flow in a given geometry, the flow is solely controlled by Re, For viscoelastic fluid flow, however, it is not so straightforward. The elastic effects are introduced via the extra stress terms, and the relative Corrections caused by elasticity and inertia depend upon the relative strengths of We and Re. If we introduce R and λ as the characteristic length and time, then the velocity and stress are scaled by R/λ and η0 /λ. In this case, the momentum equations in dimensionless form have the form: u ∂u uu − El∇u u) = El{−∇pII + ∇ · [τ − 2(1 − E)βD D ]}, + ∇ · (uu ∂t where We η0 λ El = . = Re ρR 2

(15)

(16)

From this, it can be seen that the flow is governed by the so-called elasticity number El which is a function of the material properties of the fluid and the geometry character. The flow will be dominated by the elastic effects when El ≥ 1, otherwise, it will be dominated by inertia. The initial and boundary conditions relevant to the problem can be well defined. That is, prior to time t = 0, the fluid is at rest, and at t = 0, the lid starts abruptly to rotate with constant angular velocity . At all times, no-slip condition is assumed around the confined cylinder walls and lid. As a result, meridional and azimuthal circulations develop. All of the results will be presented in dimensionless form with the length, velocity, time as well as stress being scaled by R, R,1/ as well as η0 , respectively. 3. Numerical method The initial-boundary-value confined flow problem is solved by using the efficient implicit FVM. The detailed formulation procedure of the method for Cartesian coordinates has been well documented for heat transfer and Newtonian fluid flow problems by Patankar [35], and by Xue et al. [28] for viscoelastic flow problems. In this section, we will briefly introduce the formulation procedure in a general orthogonal coordinate system, and discuss some specific difficulties arising from coordinate transformation in the solution method. 3.1. Finite volume formulation The key concept utilised throughout the finite volume formulation is the principle of conservation of a certain physical quantity expressed by the governing equations over any finite volume, or control volume,

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

345

Table 1 The definition of the coefficients and the source terms for governing equations Equation

8

3

0

S8

Continuity

1

ρ

0

0

ρ

η0

λ

0

Momentum

Constitutive

  vr u = vθ  vz  τrr     τrθ    τrz τ = τθ θ     τ    θz τzz

  Sr S 8 = Sθ  Sz  Srr     Srθ    Srz S8 = Sθ θ     S    θz Szz

or control volume. The discretisation process becomes more convenient by recognizing the fact that all of the relevant governing equations possess a common form, that is, the form of the general transport equation   1 ∂ h1 h2 h3 ∂ 3uk 8 unsteady term 38 + convection term ∂t h1 h2 h3 ∂xk hk   ∂ h1 h2 h3 ∂8 0 (17) = diffusion term + source term SΦ , ∂xk ∂xk h2k where Cartesian-tensor notation has been used. h1 = 1, h2 = r, h3 = 1, x1 = r, x2 = θ, x3 = z, and u1 = vr , u2 = vθ , u3 = vz in a cylindrical coordinate system, and h1 = h2 = h3 = 1, x1 = x, x2 = y, x3 = z, and u1 = u, u2 = v, u3 = w Cartesian coordinates. The dependent variable 8 can be a component of a vector or a tensor and even a constant. The coefficients 3, 0 have different meanings for each different dependent variable, and S8 is called the source term. which lumps all the terms that cannot be accommodated in the convection and diffusion terms, and is specific to a particular 8. The contents of the coefficients and the source terms for different dependent variables, thus, different equations governing viscoelastic fluid flows are listed in Table 1 with the details of the source terms being given in Appendix A. With the finite volume formulation, the flow domain is discretised into a set of non-overlapping control volumes, which can be irregular in size and shape (we are only concerned with structuredorthogonal-control volumes here). A section of such a discretised domain in 2D form is shown in Fig. 2a with a typical 3D control volume for node P showing in Fig. 2b. The dependent variable 8 is calculated and stored at the centroid (say, node P) of the control volumes. The algebraic equation for each variable at its node is obtained by integrating the general transport equation over the time interval δt and the respective control volume with the geometric volume 1V = h1 h2 h3 dx1 dx2 dx3 and using the divergence theorem whenever possible to reduce the dimensionality of integrations, and then introducing the appropriate temporal and spatial approximations to replace the integration:

346

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

Fig. 2. (a) Typical 3D control volume for node P, and (b) partial view of a 2D staggered mesh system.

Z Z

∂8 3 dV dt + ∂t δt 1V Z Z S8 dV dt. = δt

Z Z δt

1V

∂ h1 h2 h3 ∂xk



h1 h2 h3 ∂8 h1 h2 h3 3uk 8 − 0 hk ∂xk h2k

 dV dt (18)

1V

Being written in the form shown above, except for the source term, the governing equation is identical to the equation under Cartesian coordinates. Therefore, a complete derivation of the algebraic equation for the dependent variable 8. which has been addressed in our previous work [28], requires no further elucidation here, and the final discretised equations can be obtained by an analogous derivation. For the sake of generality, the final discretised equation, namely, an algebraic approximation equation which establishes the relation of the value 8P at the centered node P to its neighboring nodal values can

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

be symbolically expressed in a generalized form for each of the control volumes in the domain X aP 8P = anb 8nb + S¯c + aP0 8nP ,

347

(19)

nb

where the superscript n denotes the values taken from previous time level, and X 3 anb + aP0 − S¯P , aP0 = 1V ; aP = δt nb

(20)

with the summation being to be taken over all of the neighboring nodes nb of the centered node P. As usual, the volume integral of the source term S8 has been linearised as a function of the 8P with S¯P being the part that does not explicitly depend on 8P , and S¯P , the coefficient of 8P . An overbar means the applied values are evaluated using the known field at time level n (or from the previous iteration for steady-state calculations). In view of its simplicity and unconditional stability for transient calculations, the first-order backward Euler implicit formula is adopted for the temporal approximation. The coefficients anb are the functions of the dependent variables, and their contents depend on both the form of the control volume used and the approximation scheme adopted. It is these coefficients that determine the temporal and spatial accuracy of the final solution. Our previous experience [32] indicated that the power-law (PL) scheme of Patankar [35] works well for the discretisation of the convective-diffusive (momentum) equation because of the simplicity of its implementation and low computational expense as well as better conservational properties. The validity of the method is of course not restricted to this choice. With the adopted scheme and mesh system shown in Fig. 2, the coefficients anb take the form anb = Dnb f (|Pnb |) + dsign(nb)Fnb , 0e,

(21)

where the ‘sign(nb)’ is + for the upstream faces (w, s, d) and − for the downstream faces (e, n, u). The function f(|Pnb |) has the form f (|Pnb |) = d0, (1 − 0.1|Pnb |)5 e, where Fnb = (Λuk A)nb ;

 Dnb =

0A δxk

(22)  (23) nb

can be thought of as being the strength of the convection (or flow) through the control volume surface nb with area (A)nb ; and the diffusion conductance, respectively. The ratio of the two strengths is called the local Peclet number:     F 3uk δxk = , (24) Pnb = D nb 0 nb with (δxk )nb the distance between the central node P and its neighboring node in the k direction. In cylindrical coordinate system, for example, we have:    (δr)n and (δr)s in the r direction in the θ direction (δxk )nb = rm δθ   (δz) and (δz) in the z direction u d

(25)

348

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

and (A)nb

  rn 1θ 1z and rs 1θ 1z = 1r 1z  rm 1θ 1r

for surface n and s for surface e and w for surface u and d.

(26)

With the discretised form given in Eq. (19) for the general transport equation (Eq. (17)), the discretised form for a specific equation can be readily written out by imposing different contents on the coefficients 3, 0 and the discretised source terms. As far as the discretisation of the constitutive equations is concerned, since the decoupled technique is used to treat the interaction of kinematics and stresses in this work, which means that the non-Newtonian stress contribution in momentum equations is treated as a body force and the extra stresses can be calculated separately via the adopted constitutive equation with known kinematics from the previous iteration (or time level), some specific techniques can be introduced to take its hyperbolic nature into account, such as streamline upwind (SU) of Brooks amid Hugiles [36], or the strain tensor tracking algorithm developed by Luo and Tanner [37] for integral viscoelastic models [38]. However, as shown in our previous work, the finite volume formulation could be directly extended to its discretisation without any extra effort and the use of the first-order upwind difference (UD) for spatial discretisation effectively ensure numerical stability of the resulting discretised system. The discretised equation obtained in this manner expresses the conservation principle for 8 in the finite volume. The most attractive feature of the control-volume formulation is that the resulting solution would imply that the integral conservation of quantities such as mass, momentum, and energy is exactly satisfied over any group of control volumes and, of course over the whole computational domain. This characteristic exists for any number of grid points, not just in a limiting sense. Thus, even the coarse-grid solution can exhibit the exact integral balances. 3.2. Solution method After implementing the above mentioned steps for each of the control volumes in the domain, then each of the nodes in the flow domain has its discretised equation like Eq. (19). The final set of nominally linear algebraic equations for each of the nodes in a 3D domain include three momentum equations for the velocity vector uk (k = 1, 2, 3), six constitutive equations for the stress tensor τ ij as well as the continuity constraint. The non-linear coupling, such as the coupling of the kinematics and dynamics (stress), implicitly appears through the coefficients and source terms. With this primitive variable formulation, the equation for evaluating the pressure in which the velocity and the pressure are coupled can be provided by enforcing mass continuity over each cell, which has been described elsewhere and therefore does not require further explanation here. Without a special reconstruction of the control volume surface velocities (see, e.g. Rhie and Chow [39]) to prevent spurious oscillations in the pressure field, a staggered mesh system is used, that is, the control volume for the velocity component in the specific direction lies between the pressure points in the respective direction, as shown in Fig. 2b. The control volume for other scalar or tensor variables is the same as that for the pressure. It seems that the mesh system and the form of the algebraic equation for the dependent variables are completely analogous to those in Cartesian coordinate system. Thus, the solution methodology used there may directly be employed here.

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

349

However, with the cylindrical coordinate system (or spherical coordinates), there are some numerical difficulties related with the specification of boundary conditions. especially for the calculations with a structured mesh system with the axis or centre being included in the solution domain as a boundary and without the assumption of axisymmetry. One of the difficulties is associated with the geometry singularity on the axis where r = 0 and cells around the axis degenerate from quadrilateral ones to triangular ones. Terms containing 1/r tend to become singular and terms like ∂8/∂ ␪ become undefined. Therefore, in order to solve a non-axisymmetric flow problem involving a hollow cylinder or sphere, special treatment is required on the axis. To eliminate this singularity, various methods have been proposed, such as the use of a mixed Cartesian (for the region near the axis) and cylindrical (for the region away from the axis) grid system, or use of a change in variable [40], or simply omitting a mesh point at the axis from the set of mesh points [41]. Recently, in developing a FVM in polar cylindrical non-staggered grid, Inghami and Wen [42] employed a circular cell on the axis to replace the triangular cells around the axis. However, although much extra effort has been devoted in discretisation, the numerical results seem not very satisfactory even for axisymmetric flow problems. Some simple, but effective technique to circumvent this singularity is proposed by Ozoe and Toh [43]. They simply interpolated a centre radial velocity component ur (0,θ) by using the two velocity vectors, v (1r, θ ± 0.5π ) on its adjacent grid points, where 1r is the radial distance of the circle or cylinder (in 3D) composed of the triangular cells around the axis, that is, ur (0, θ ) =

vθ (1r, θ + 0.5π ) + vθ (1r, θ − 0.5π) . 2

(27)

This simple procedure is followed in this work. Actually, if the axis is taken as one of the radial boundaries, what we need are the centre values of the radial velocity component. Once these values are determined, the calculation can be carried out, and the centre values of the vertical axis velocity component w(0,θ,z) (for 3D), as well as the stress tensor or any scalar variable can be directly interpolated approximately using summation over surrounding points. Thus, no calculation is needed for those triangular cells around the centre and no terms containing 1/r, or ∂8/∂θ are needed in the w hole calculations. With the cylindrical coordinate system, another difficulty is that there is no definite boundary in the circumferential direction. In our work, this is effectively circumvented by using the cyclic tridiagonal matrix algorithm which is detailed in the Appendix B. After removing the above mentioned obstructions, the solution methodology for such a set of discretised equations involving multi-dimensions and multi-couplings developed in Cartesian coordinate system can be directly employed here. In this work, we employ the SIMPLEST (SIMPLE [35] with Splitting Technique) [44] developed in our previous work by combining the advantages of the SIMPLER (SIMPLE Revised) [35] and the PISO (Pressure-Implicit with Splitting of Operators) schemes [45].

4. Numerical results Before carrying out the simulations aiming at reproducing the experimental observations, for the steady-state Newtonian solution, the accuracy of the code developed is checked by comparing the predicted results with the experimental point-wise measurements with various grid sizes on vertical diametral

350

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

plane. Based on these comparisons, the mesh resolution required to obtain a size-independent solution is determined for the rest of simulations. Then a time-step employed to predict the transient process is determined by checking the repeatability of the solutions on a given resolution for more than two different time-steps. As far as the convergence criterion for calculations is concerned, the resulting field will be regraded as steady flow for transient calculations and as convergent solutions for steady-state calculations when the integral dimensionless residual for each of the governing equations over all of the control volumes tend to zero (about O(10−10 ), machine accuracy) with further increasing times/number of iterations. Overall, the simulations will concentrate on the following four points: 1. Numerical verification of part of Escudier’s visualisation [7] on Newtonian fluids, thus checking the performance of the code further. 2. Numerically checking the observations of Spohn et al. [14] about the opened and asymmetric bubble structure in the steady vortex breakdown process. 3. Numerical verification of the observed conclusions of Stokes et al. [20] on the shift of the existance domain of vortex breakdown for a dilute flexible polymer solution by using the UCM model, thus evaluating the suitability of the model in the flow process.. 4. Numerical investigation on the double-cell flow structure observed by Escudier and Cullen [17] for a strong shear thinning viscoelastic liquid by using the SPTT model, thus checking model’s predictability.

4.1. Steady-state and transient Newtonian flows Experimentally, using an LDA technique, Michelsen [46] measured the point-wise distributions of velocity components vθ and vr along the lines with constant r in a confined cylinder with H/R = 1.0 rotating at a Reynolds number of 1800. In Figs. 3 and 4, numerical predications with three mesh resolutions (Ml: 1r × 1z = 1/20, M2: 1/45 and M3: 1/60) on a vertical diametral (rz) plane are compared with the measurements at two of the positions, r = 0.6 and r = 0.9, respectively. As we can see, M3 is fine enough to reproduce the measurements. With M3, one of the typical visualisations observed by Escudier [7] for a Newtonian fluid in a confined cylinder with H/R = 2 at Re = 1854 is numerically reproduced as shown in Fig. 1b. As we can see, the predicted flow structure is perfectly axisymmetric. To view the results more clearly, in the rest of the figures showing the vortex breakdown structure only half of the 2D vertical diametral plane with the left vertical boundary being the axis of symmetry will be presented. In Fig. 5a–c, three of typical visualisation observed by Escudier [7] for a Newtonian fluid in a confined cylinder with H/R = 2 at Re = 1002, 1492 and 1854 are presented. As seen, both the entire flow structure and the bubble positions have been predicted quantitatively. With time-dependent calculation, it is possible to capture numerically the transient flow process of the development of vortex breakdown. With our implicit temporal discretisation method, there is no obvious time-step restrictions, such as the CFL condition for explicit formulation. However, there does exist an optimal time-step at which the real physical history of flow can be captured efficiently. There is no set rule to obtain this optimal value, which may depends on the mesh resolution used and value of Re as well as value of We, for viscoelastic calculations. Usually one has to do some numerical experiments to check the repeatability of the solutions obtained with various time-step in a given mesh solutions.

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

351

Fig. 3. Profiles of rotational velocity vθ and radial velocity vr along r = 0.6.

However, some light can be shed on finding out the time-step range for efficient transient flow calculations. As we know, the steady-state solution can be reached by using the time-marching scheme in which a positive inertia is introduced via the terms related with the time-step, i.e., ap0 in Eqs. (19) and (20), and it will play as a local under-relaxation function. In this case, the under-relaxation factor 2 in the case of S¯p = 0 is P anb 1 = , (28) 2 = P nb 0 1+κ nb anb + ap with, κ=P

ap0

nb

anb

From this relation, by assuming a uniform dimensionless mesh size 1h, and the uniform material properties 3, 0 in all directions, for the case when the flow is dominated by convection as for the swirling flow, the dimensionless time step δt would roughly be δt ∼

1h . κ(u1 + u2 + u3 )

(29)

Since in the swirling flow, the dominant velocity is u2 (=uθ = 1) and u1 (=ur ), u3 (=uz )  1, we have δt ≤

2 1h 1h = . κ 1−2

(30)

352

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

Fig. 4. Profiles of rotational velocity vθ and radial velocity vr along r = 0.9.

Fig. 5. Numerically simulated streamlines at steady state in a confined cylinder with H/R = 2 and (a) Re = 1002, (b) 1492, and (c) 1854 for Newtonian fluid.

For steady-state calculations, for a given mesh resolution, the value of 2 depends on Re for Newtonian flow, as well as We for viscoelastic flow. For example, for the above steady-state calculations, the values of Re range from 1002 to 1854, the corresponding values of 2 are in 0.85–0.75 if M3 is used. Thus, for the corresponding transient calculations, δt will be roughly in the range of (5.67–3)1h (0.0945–0.05 dimensionless units). This range is by no means the limitation on the time-step for transient calculations. In the above analysis, some stricter assumptions, such as a uniform mesh in terms of the minimum size in all directions and S¯p = 0, have been made. Thus, the actual time-step for transient calculations

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

353

Fig. 6. Evolution of flow field in time for Newtonian fluid in a confined cylinder with H/R = 2 and Re = 1854 predicted using (a) δt = 2.76 1h, and (b) δt = 22.1 1h on a given mesh M3.

immune from oscillations leading to a unphysical solution should be larger than the upper limit of this range. We select the case (Re, H/R) = (1854, 2) to carry out Newtonian transient calculations to demonstrate numerically the development of flow field at different times, thus the evolution process of vortex breakdown in time, after the start of the boundary rotation. First. we use a time step smaller than the upper limit of the above range, δ = 2.76 1h, in the simulation. The predicted evolution of flow field in time is shown in Fig. 6a. As seen, the swirling flow field experiences a complex flow process before it reaches the steady state. Before t = 90, the streamlines are seen to be monotonic with the swirling centre moving down gradually to the rotating edge. At about t = 110, a bubble first appears on the axis of symmetry. This bubble quickly changes shape and moves down, fol-

354

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

Fig. 7. Evolution of the axial velocity component vz in time at a point z = 1 along the centre line of the confined cylinder.

lowed by another small bubble appearing at t = 130. Then, at t = 150, just one large bubble is left. However, at t = 190, this bubble is split into two small bubbles. Finally, when the steady state is reached at t = 350, two bubbles move closer with the upper one growing and the lower one diminishing. This simulation is repeated with another quite large time step, that is, δt = 22.l1h, as shown in Fig. 6b. The agreement between the two set of predictions implies that the time step 22.l1h is small enough to capture the real physical history of flow for this case, and with the implicit formulation, the minimum time step required for convergent solution is quite large compared to that required in an explicit formulation. Sorensen and Loc [10] simulated the same problem by using an explicit finite difference method with a combination of fourth-order and second-order accurate schemes based on a streamfunction-vorticity-circulation formulation for axisymmetric flows. They failed to capture a convergent solution after δt > l51h. Judging from the predicted position of the stagnation point of the large bubble in the steady state, we can see that our results are closer to the experimental visualisation, while Sorensen and Loc’s predictions differ from the experimental visualisation by about 3% of the cylinder height. That the solution is independent of the time step after the time step is reduced to some level is further demonstrated in Fig. 7. where the evolution of the axial velocity component in time at a point z = 1 along the centre line of the confined cylinder is compared for different time-step. As seen, for the cases when δt ≤ 28.21h, all of the solutions are almost overlapped, and at about δt = 30.01h, a unphysical solution history is produced. Also, from this figure, we can see that a steady flow regime is established at about t = 300–350. Based on data files obtained from the 3D numerical simulations without the axisymmetric assumption, it is convenient to investigate the flow features, such as the structure of the bubbles in vortex breakdown, by tracing the flow path of a fluid element in 3D space. One of the interesting observations recently made by Spohn et al. [14] is that even in steady flow, the flow is not axisymmetric in such an axisymmetric

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

355

Fig. 8. Numerically predicted trace of fluid elements around the breakdown bubbles in terms of 3D streamlines for a Newtonian fluid in a confined cylinder with H/R = 1.75 and Re = 1850.

geometry. Numerically, we mimic their experiment by choosing the fluid elements symmetrically located around the centre axis near the upper rigid cover, and tracing their paths in the confined swirling flow. In Fig. 8a, four fluid elements along the line A–O–B are symmetrically chosen. As seen, the symmetrical elements flow dowmstream helically and symmetrically. Thus is further demonstrated in Fig. 8c by viewing from the top of the geometry. From such a flow structure, we also can see that the breakdown bubbles in 3D space are formed by a series of extension spiral spring-like streamlines, thus they are open with inflow and outflow as observed by Spohn et al. [14], but they are highly axisymmetric from the upstream side of the bubble to their real side, which is different from that observed. For this discrepancy, one of the explanations may be the small imperfections of the experimental rig. Another possibility is that although the kinematic fields are axisymmetric, the integral curves derived from the kinematic fields may not be axisymmetric. This needs further investigations as most flow visualisation represent integral curves. To show how complex the swirling flow in 3D space is, we trace the path of a fluid element near the edge of the rotating lid and present it in Fig. 8b. From this, we can see that the vortices on the vertical diameteral plane, as shown in Fig. 8d, are actually composed by a series of spiral spring-like, swirling streamlines. 4.2. UCM fluid The above numerical results have demonstrated that the confined swirling flow of Newtonian fluid in a given cylinder is controlled solely by inertia. Now we start investigating the purely elastic effects on the swirling flow by including the non-linear constitutive equations used to describe viscoelastic fluids with a constant viscosity in calculations. From the above dimensional analysis, we have seen that for viscoelastic flows, if El  1. then inertia still dominates the flow field, that is, a steady Newtonian-like flow is expected. However, the flow field will

356

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

Table 2 Vortex breakdown feature for different fluids Cylinder H/R

Reynolds Number (Re)

Number of breakdowna Newtonian

25 ppm

45 ppm

1.5

1139 2000

◦ ×

◦ ×

× ◦

2.0

1500 2080

◦ ◦◦

◦ ◦◦

× ◦

2.5

1918 2494

◦ ◦◦

× ◦◦

× ◦◦

a

One bubble — ◦, two bubbles — ◦◦ and no bubble —×.

be dramatically influenced due to elasticity, which has been demonstrated in the series of experimental work by Stokes et al. [19–21] using polyacrylamide (PAA) Boger fluids with a constant viscosity. One of the conclusions made by observing the swirling flow of low viscosity dilute PAA solutions is that for a given H/R, a larger values of Re for the appearance and disappearance of vortex breakdown bubbles is required compared to a Newtonian fluid of similar shear viscosity. The values of Re, required will increase with increasing polymer concentrations, and to some extent vortex breakdown appeared to be totally suppressed at the maximum Re tested. Numerically, we will employ the UCM model to describe the low viscosity dilute PAA solutions used in the experiment, and the polymer concentrations will be characterized by the relaxation time Although the exact data of λ for the solutions are not available at the time, roughly, they are about 0.007 and 0.026 s for 25 and 45 ppm PAA solutions with the shear viscosity being about 41 and 46 mPa s, respectively [47]. To numerically verify the experimental observations, a series of numerical experiments have been carried out using the UCM model with the material parameters being η0 = 41 mPa s and λ = 0.007 and 0.020 s for 25 and 45 ppm PAA solutions, respectively. The typical results are listed in Table 2 where the appearance and disappearance of vortex breakdown bubbles in three different cylinders with H/R being 1.5, 2.0 and 2.5 at Re ranging from 1139 to 2494 for three kinds of fluids, Newtonian, 25 ppm and 45 ppm PAA solutions, are tabulated. Similar to the experimental observations, the onset of vortex breakdown and its disappearance are delayed due to elastic effects, and the higher the concentration, the more remarkable the delay. For example, in a cylinder with H/R = 1.5, at Re = 2000, there is no breakdown bubble for Newtonian fluid, as well as for 25 ppm PAA solution. However, one breakdown bubble is in evidence for the PAA solution with higher polymer concentration (45 ppm PAA). Such a delay is further demonstrated in Figs. 9 and 10. As shown in Fig. 9a–c in cylinder H/R = 1.5, at Re = 1139, one breakdown bubble appears for Newtonian fluid. This breakdown bubble is obviously weakened with its position moving down for 25 ppm PAA solutions, and it totally disappears for 45 ppm PAA solution. In Fig. 10a–c, the vortex breakdown bubbles for 25 ppm, and 45 ppm PAA solutions as well as the corresponding Newtonian fluid in H/R = 2.5 at Re = 2494 are presented. As seen, for the Newtonian case the visualisation observation presented by Escudier [7] in their Fig. 4e has been numerically reproduced here, while the flow pattern at the same conditions for 25 ppm PAA solution described by the UCM model does not seem to be different. For our modelled 45 ppm PAA solution, a smaller value of λ (=0.02 s) than reported (0.026 s) has been used. As a result, the predicted flow structure is not in the one breakdown bubble region as depicted in the

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

357

Fig. 9. Numerically predicted vortex breakdown bubbles for (a) 25 ppm PAA solution and (b) 45 ppm PAA solution, as well as (c) the corresponding Newtonian fluid at H/R = 1.5, Re = 1139.

Fig. 10. Numerically predicted vortex breakdown bubbles for (a) 25 ppm PAA solution and (b) 45 ppm PAA solution, as well as (c) the corresponding Newtonian fluid at H/R = 2.5, Re = 2494.

experimental work of Stokes et al. [20] Fig. 2. Instead, there are still two bubbles. However, the flow structure is more or less like the case of Newtonian at a lower Re (Re = 2126) shown in Escudier’s Fig. 4d in [7] in which the positions of the breakdown bubbles are lower, the sizes are smaller, and the distance between two bubbles are larger. All of those indicate that the vortex breakdown process is dramatically delayed in terms of Re for a given geometry. Obviously. the slight elastic effects on Newtonian flow field for viscoelastic fluids with a low viscosity are due to the introduction of the first normal stress difference N1 which can be characterized by the elastic number El. It is N1 that counterbalances the inertia force, thus delaying and even suppressing vortex breakdown. For viscoelastic fluids with a low viscosity, El is quite small (in the range of 0.01–0.1 in above calculations), so the flow field is Newtonian-like. However, different flow patterns are expected with increasing values of El due to the stronger elastic effects. The specific pattern for a given fluid will depend on the balance between the elastic force (El or N1 ) and the inertia force (Re). Once N1 is large enough to counterbalance the inertia force, then the

358

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

Fig. 11. Dramatic difference in steady swirling flow patterns in terms of streamlines in the same confined cylinder with H/R = 1 at Re = 0.32 for: (a) a Newtonian fluid; and Boger fluids at (b) λ = 0.05 (We = 0.006), (c) λ = 0.1 (We = 0.01), (d) λ = 0.15 (We = 0.02), (e) λ = 0.36 (We = 0.05), and in the low row, the corresponding close-up view of the meshes and flow field in terms of streamlines around the edge of the rigid cover.

flow will be dominated by the elastic normal stress, and the flow field will be partially or even entirely reversed depending on the values of El. Stokes et al. [19] have shown the totally reversed steady flow in their visualisation for a Boger fluid with λ = 0.36 (We = 0.05) at Re = 0.32 (thus, El = 0.15) in H/R = 1. As shown in Fig. 11e, the same flow pattern has been successfully predicted numerically with the UCM model under the experimental conditions. Our numerical experiments at the same value of Re, but different values of λ, thus, different elastic fluids, demonstrated that at this low value of Re, the totally reversed flow would not occur until the elasticity of the fluid reaches to some high level. As shown in Fig. 11a–c, the flow patterns for fluids with different λ (from 0 to 0.36. thus El = 0 to 0.15) are presented. We can clearly see the competitive process between the elastic force and the inertia force. A reversed elastic ‘ring’ vortex is developed at the edge of the rotating lid for the fluid with λ being about 0.05, and with increasing values of λ, this ring vortex is growing and pushing Newtonian vortex to become a central ‘ring’ vortex, and with further increase of λ (about 0.2), this central ring vortex is diminished and the entire field is occupied by one reversed elastic vortex. This process has actually been observed quantitatively in the experiments by Hill [15]. A similar process is also predicted in the geometry with H/R = 2.0 as shown in Fig. 12a–d. We notice that the only difference is that there is a small and weak ring vortex appears at the edge of the top rigid cover of the cylinder at this low value of Re. A close-up view of the structure of this ring vortex is shown in the lower row of the same figure, we find that it occupies several meshes and is counter-rotating to the neighboring larger vortex. In our calculations, we always found such an edge ring vortex when the flow field is not dominated by inertia (at a low or medium high value of Re) in a geometry when H/R ≥ 2 and its size could increase with increasing H/R. However, it disappears when inertia dominates the flow. Due to being small in size and weak in strength of the edge ring vortex, it may be difficult to be observed in experimental visualizations. Overall, our numerical experiments indicate that the experimental observations on the swirling flow can be reproduced numerically, which clearly indicates that the code is working well for viscoelastic flow calculations.

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

359

Fig. 12. Dramatic difference in steady swirling flow patterns in terms of streamlines in the same confirmed cylinder with H/R = 2 at Re = 0.32 for: (a) Newtonian fluid; and Boger fluid at (b) λ = 0.05 (We = 0.006), (c) λ = 0.1 (We = 0.01), (d) λ = 0.15 (We = 0.02), (e) λ = 0.36 (We = 0.05). Table 3 Model parameters for the test fluid Parameters

η0 (Pa s)

β

ρ (kg/m3 )

␭ (s)

ε

Values

1.91

0.99

1200

0.016

0.12

4.3. PTT fluid In this section, we will numerically investigate another one of the interesting phenomena observed by Escudier and Cullen [17] in the confined swirling flow for viscoelastic fluids with a strong shear thinning viscosity, but low elasticity compared to the inertia force (El  1), that is, a double-cell flow structure. In experiments. aqueous solutions of polymer carboxymethylcellulose (CMC) with concentrations ranging from 0.75% to 1.5% w/w were used, and the strong shear thinning viscosity was well represented by the Cross model. We simulated one of the experimental cases shown in the Fig. 5 with (H/R, Re, El) = (2.0, 67.8. 0.0029). The test fluid is modelled by the single-mode SPTT constitutive equation with the numerical parameters listed in Table 3 based on their Table 1 [17]. Here, we must note that the Weissenberg number (We) defined in their work is the elastic number (El) used in this work. In our simulations, no attempt has been made to match exactly the shear thinning viscosity. Instead, the parameters  = 0.12 and ␭ = 0.0016 are used to characterize fluid shear-thinning and elastic behavior in the SPTT model. The relaxation value λ is worked out based on the relation of El and λ, that is, λ=

η0 El . ρR 2

(31)

We also note that a shear dependent elastic number, in which the relaxation time is defined as N1 /2η0 γ as shown in their Fig. 2, was used to characterize the flow conditions in their experiments as shown in

360

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

Fig. 13. Steady swirling flow patterns in terms of streamlines in a confined cylinder with H/R = 2 for strong shear-thinning viscoelastic fluids described by the SPTT model at Re = 67.8, with (a) λ = 0 (Newtonian), (b) λ = 0.016 (El = 0.0029), (c) λ = 0.08 (El = 0.014), (d) λ = 0.15 (El = 0.027), (e) λ = 0.2 (El = 0.036), (f) λ = 0.25 (El = 0.045), (g) λ = 0.3 (El = 0.054), (h) λ = 0.4 (El = 0.072).

their Fig. 3. However, in numerical calculations, with the constitutive model, a zero-shear rate relaxation time is usually used. Therefore, it should be much higher for strong shear thinning viscoelastic fluids. Our simulate results indicate that with the above fit, the double-cell structure cannot be predicted. In Fig. 13a, the corresponding Newtonian flow field is presented. Again, at this medium high Re at which the flow field is neither dominated by inertia or elasticity, and no vortex breakdown occurs (roughly, Re is of order O(1–500)), the edge ring vortex as discussed in the previous section is in evidence. At the prescribed value of ␭ (=0.016 s) (see, Fig. 13b), the flow is still Newtonian-like with a small reversed edge ring vortex there, in contrast to the observed double-cell structure. Near the rotating lid, shear rates reach the maximum, the swirling velocity is higher due to the strong shear thinning effects, thus the flow in this region is dominated by the inertia force (Newtonian-like). In the region far away from the rotating lid, the shear rates are low in view of the weak flow there, thus the flow is dominated by the elastic effects (reversed flow) even for the fluid with low elasticity. Therefore, a double-cell flow structure should be expected. With the SPTT model, a zero second normal stress difference is predicted. The role of the second normal stress difference was briefly investigated with the PTT model with ξ = 0.1 (N2 /N1 = −ξ /2), as well as the Giesekus model [48] with the mobility factor α = 0.1. However, the similar results were obtained. We agree with the conclusion made by Kramer [22] that the influence of the second normal stress difference is not of significance. Then, the only possibility is the discrepancy in the value of λ which controls the intensity of the first normal stress difference. As noted above, in their experiments [17], a shear dependent relaxation time was used to difine their elastic number. If a zero-shear rate relaxation time was used, the elastic number would be about 30 times higher according to their Fig. 4 when Re → 0. Therefore, the zero-shear rate relaxation time would be in the range of about 0.25–0.3 s for the test fluid according to their Fig. 2. We continue our simulations by gradually increasing values of λ. As shown in Fig. 13c–h, the

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

361

reversed ring vortex starts growing with incrasing values of λ, and at. λ = 0.25, the observed double-cell flow pattern is in evidence. With further increasing values of λ, the reversed vortex occupies more and more space and the Newtonian vortex is compressed to the edge of the rotating lid. And, finally, the entire flow field may reversed or become unstable. However, we stop our calculations at λ = 0.4. The value of We at λ = 0.4 is about 4.8 or elastic number is about El = 0.07 which is about 25 times higher than that used in the experiment when a shear dependent relaxation time was used. Finally, it is worth mentioning another interesting issue. It is noted that in most of experimental work, a central ring vortex around the centre of the rotating lid was reported for both Boger fluids (see, e.g., [21]), and shear thinning viscoelastic fluids (see, e.g., [16]) before the flow field becomes unstable in a given geometry. By closely examining the simulated results shown in Fig. 13, we note that when λ ≥ 0.25, a small central ring vortex around the centre of the rotating lid is in evidence. This may be a precursor of the impending instability. We will continue our numerical investigation on this issue in the next part of the work with time dependent calculations, and hopefully, we can numerically check the hypothesis about the stability criterion (see, e.g. [1]).

5. Conclusions In this work, aiming at numerically reproducing the steady and/or unsteady complex Newtonian or viscoelastie flows observed in a confined cylinder with a rotating lid, an efficient implicit FVM in a general orthogonal coordinate system, particularly, for a cylindrical coordinate system without the assumption of axisymmetry has been developed. For steady Newtonian flows, excellent quantitative agreement between predictions and experimental observation has been reached for such convective-dominated flows (high Re number), thus, confirming the accuracy and efficiency of the numerical method. For the implicit temporal discretisation method, some light has been shed on the optimal time-step range for stable transient flow calculations. By carrying out time-dependent calculations to simulate the transient flow process of development of vortex breakdown in the flow domain, it is found that the minimum time step required for convergent solution can be quite large compared to that required in an explicit formulation. The flow patterns of vortex breakdown in 3D space have been clearly presented based on data files obtained with the 3D numerical simulations without the axisymmetric assumption. It is found that the bubbles around the centre axis are formed by a series of extension spiral spring-like streamlines, thus they are open with inflow and outflow as observed in experiments, but they are highly axisymmetric from the upstream side of the bubble to their rear side, which is different from some observations. The numerical results with the UCM model confirmed the experimental observations for some Boger fluids, that is, the onset of vortex breakdown and its disappearance will be delayed due to the slight elastic effects, and the higher the concentration, the more remarkable the delay. Different flow patterns are predicted with increasing values of El. Depending on the values of El, the flow field will be partially or even entirely reversed, and finally become unsteady, periodic, and non-axisymmetric. The steady double-cell flow structure for viscoelastic fluids with a strong shear thinning viscosity, but low elasticity compared to the inertia force (El  1) has been successfully predicted with the test fluids being modeled by the SPTT.

362

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

From this study, we have demonstrated that the implicit FVM developed can be efficiently employed for time-dependent calculations. Continues work focusing on simulating unsteady viscoelastic flows will be reported in the near future. Acknowledgements This work is supported by the Australian Research Council. The calculations were done on the comupting facilities of the Sydney Distributed Computing Laboratory (SyDCom). The suppport is gratefully acknowledged. Appendix A. Contents of source terms In this appendix, the source terms for each of the governing equations with (E = 1) or without (E = 0) the EVSS formulation are listed. 1. The momentum equations for three velocity components vr , vθ and vz (for Newtonian fluids, β = 0)      vr ∂ ∂vr ∂vθ ∂p 1 2 + ρvθ − η0 + − η0 Sr = − ∂r r r∂θ r r∂θ r      ∂τθr 1 ∂ ∂τrz τθθ ∂ r∂vr + (rτrr ) + + − − (1 − E)β η0 r ∂r r∂θ ∂z r r∂r ∂r        vθ ∂ η0 ∂vθ vr ∂vr ∂vr ∂ η0 − + η0 − + (A.1) + r∂θ r∂θ r ∂z ∂z r r∂θ r    ∂p 1 vθ ∂  vr  ∂vr − ρvθ vr − η0 − + η0 Sθ = − r∂θ r r∂θ r r∂θ r   1 ∂ ∂τθθ ∂τθz τrθ + (rτrθ ) + + + r ∂r r∂θ ∂z r         r∂vθ ∂vθ ∂vθ ∂ ∂ vθ ∂ −(1 − E)β η0 + η0 + + η0 r∂r ∂r r∂θ r∂θ r ∂z ∂z   vθ η0 ∂vr − + r r∂θ r   ∂p ∂τθz ∂τzz 1 ∂ + (rτrz ) + + Sz = − ∂z r ∂r r∂θ ∂z        ∂ ∂ r∂vz ∂vz ∂vz ∂ η0 + η0 + η0 −(1 − E)β r∂r ∂r r∂θ r∂θ ∂z ∂z

(A.2)

(A.3)

2. The PTT constitutive model for six stress components τ rr , τ θθ , τ zz , τ rθ , τ θz , and τ rz (For the Oldroyd-B model, set g = 1 and ξ = 0, and further, for the UCM model, β = 1):

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367



363


∂vr ∂vr − gτrr + 2λ (1 − ξ ) τrr + EτrrN ∂r ∂r       

∂vr ξ ∂vr ∂  vθ  ∂vz ∂vr ξ ∂vr N N − +r τrθ + Erθ + − + τrz +Eτrz + r∂θ 2 r∂θ ∂r r ∂z 2 ∂z ∂r  
1 ∂


∂ ∂ ∂ N N N N (A.4) λτrr + λrvr τrr + λvθ τrr + λvz τrr −E ∂t r ∂r r∂θ ∂z

Srr = 2βη0 (1 − Eg)



   
∂vθ ∂vθ vr vr N Sθθ = 2βη0 (1 − Eg) + − gτθθ + 2λ (1 − ξ ) + τθθ + Eτθθ r∂θ r r∂θ r   
∂  vθ  ∂vr ∂  vθ  ξ N − r + τrθ + Erθ + r ∂r r 2 ∂r r r∂θ    
ξ ∂vθ ∂vz ∂vθ N − + τθz + Eτθz + ∂z 2 ∂z r∂θ  
1 ∂


∂ ∂ ∂ N N N N λτθθ + λrvr τθθ + λvθ τθθ + λvz τθθ −E ∂t r ∂r r∂θ ∂z

(A.5)


∂vz ∂vz N Szz = 2βη0 (1 − Eg) − gτzz + 2λ (1 − ξ ) τzz + Eτzz ∂z ∂Z       

∂vz ξ ∂vz ∂vθ ∂vz ξ ∂vz ∂vr N N − + τrz + Eτrz + − + τθz + Eτθ z + ∂r 2 ∂r ∂z r∂θ 2 r∂θ ∂z  
1 ∂


∂ ∂ ∂ N N N N (A.6) λτzz + λrvr τzz + λvθ τzz + λvz τzz −E ∂t r ∂r r∂θ ∂z   
∂  vθ  ∂vr ∂vz N Srθ = βη0 (1 − Eg) r + − gτrθ + λ −(1 − ξ ) τrθ + Eτrθ ∂r r r∂θ ∂z      
∂vr ∂  vθ  ∂vr ξ ∂vr ∂  vθ  ∂  vθ  ξ N − r + τrr +Eτrr + − +r + r ∂r r 2 ∂r r r∂θ r∂θ 2 r∂θ ∂r r   

∂vθ ξ ∂vθ ∂vz N N − + τrz + Eτrz × τθθ + Eτθθ + ∂z 2 ∂z r∂θ    
ξ ∂vr ∂vz ∂vr N − + τθz + Eτθz + ∂z 2 ∂z ∂r  
1 ∂


∂ ∂ ∂ N N N N (A.7) λτrθ + λrvr τrθ + λvθ τrθ + λvz τrθ −E ∂t r ∂r r∂θ ∂z

364

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367



 
∂vθ ∂vr ∂vz N Sθz = βη0 (1 − Eg) + − gτθz + λ −(1 − ξ ) τθz + Eτθz ∂z r∂θ ∂r      

∂vz ∂vθ ξ ∂vz ∂vθ ξ ∂vθ ∂vz N N + + − + τθθ + Eτθθ − + τzz + Eτzz r∂θ 2 r∂θ ∂z ∂z 2 ∂z r∂θ        
∂ vθ ξ ∂  vθ  ∂vr ∂vz ξ ∂vz ∂vr N − + τrθ + Eτrθ + r − r + + ∂r 2 ∂r ∂z ∂r r 2 ∂r r r∂θ  

1 ∂


∂ ∂ ∂ N N N N N λτrz + λrvr τrz + λvθ τrz + λvz τrz × τrz + Eτrz − E ∂t r ∂r r∂θ ∂z (A.8) 

    ∂vz ∂vr ∂vr ∂vz N ) + − gτrz + λ (1 − ξ ) + (τrz + Eτrz Srz = βη0 (1 − Eg) ∂r ∂z ∂r ∂z      

∂vz ξ ∂vz ∂vr ∂vr ξ ∂vr ∂vz N N + − + τrr + Eτrr + − + τzz + Eτzz ∂r 2 ∂r ∂z ∂z 2 ∂z ∂r      
∂vz ∂vr ξ ∂vz ∂vθ ξ ∂vr ∂  vθ  N + + − + τrθ + Eτrθ − +r r∂θ 2 r∂θ ∂z r∂θ 2 r∂θ ∂r r   

1 ∂


∂ ∂ ∂ N N N N N (A.9) λτrz + λrvr τrz + λvθ τrz + λvz τrz × τθz + Eτθz −E ∂t r ∂r r∂θ ∂z Appendix B. The cyclic tridiagonal matrix algorithm In solving the generalized form of the discretised equation (Eq. (19)) with the line-by-line technique, the equation for the each node i along a line in the specific direction takes the form: di 8i = ai 8i+1 + bi 8i−1 + ci .

(B.1)

with i = 1,2,. . . ,N, where N is the maximum node number in the direction. For this set of equations, if aN = b1 = 0, which means that there are defined boundaries on both ends, the algorithms for the solution of such system are well known, such as tridiagonal matrix algorithm (TDMA). However, a so-called cyclic tridiagonal system with aN 6= b1 6= 0 will arise if an elliptic equation is discretised over the domain with periodic boundary conditions, such as, in a cylinder in which there are no defined boundaries on the circular plane, or in another word, the ’inlet’ and ‘outlet’ are overlapped. For such a system, a special algorithm called cyclic tridiagonal matrix algorithm has been developed, which can be found in Ahlberg et al. [49]. The periodic boundary conditions are expressed as: 8 0 = 8N ,

and

8N+1 = 81 .

(B.2)

The strategy to solve the system is to employ the first N − 1 equations to express 81 ,. . . ,8N−1 in terms of 8N and then determine 8N from the last equation. Suppose, we have 8i−1 = Ei−1 8i + Fi−1 8N + Gi−1 ,

(B.3)

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

365

We try to find out the relation 8i = Ei 8i+1 + Fi 8N + Gi ,

(B.4)

Where E1 = a1 /d1 , F1 = b1 /d1 , and G1 = c1 /d1 , For i = 2,3,. . . ,N − 1. By substituting Eq. (B.3) into Eq. (B.1) and rearranging, it is found out that the coefficients of equation Eq. (B.4) have the forms: Ei =

ai , di − bi Ei−1

Fi =

bi Fi−1 , di − bi Ei−1

Gi =

ci + bi Gi−1 . di − bi Ei−1

(B.5)

Then, we try to find the value of 8N . By writing Eq. (B.1) for i = N, and using condition (B.2), we have dN 8N = aN 81 + bN 8N−1 + cN .

(B.6)

Now, we try to find 81 in terms of 8N . From the first N = 1 equations, we suppose: 8i = Ti 8N + Si

(i = 1, 2, . . . , N − 1).

(B.7)

By substituting it into Eq. (B.4), we could get Ti = Ei Ti+1 + Fi

and

Si+1 = Ei Si + Gi

(B.8)

Obviously, TN = 1 and SN = 0. Thus, we can determine TN−1 ,. . . ,T1 and SN−1 ,. . . ,S1 . Finally, evaluate 8N from the Eq. (B.6) 8N =

aN S1 + bN SN−1 + cN . dN − aN T1 − bN TN−1

(B.9)

Then, the values of 8i for i = N − 1,. . . ,1 can be obtained by back substitution using Eq. (B.7). References [1] G.H. Mckinley, P. Pakdeland, A. Öztekin, Rheological and geometric scaling of purely elastic flow instabilities, J. Non-Newtonian Fluid Mech. 67 (1996) 19–47. [2] R. Keunings, On the high Weissenberg number problem, J. Non-Newtonian Fluid Mech. 20 (1986) 209–226. [3] S.-C. Xue, R.I. Tanner, N. Phan-Thien, Three-dimensional numerical simulations of viscoelastic flows-predictability and accuracy, Comput. Meth. Appl. Mech. Eng., 1999, in press. [4] H.U. Vogel, Experimentelle Ergebnisse über die laminare Strömung in einem zylindrischen Gehäuse mit darin rotierender Scheibe, Max-Planck-Institut für Strömungsforschung, Göttingen, Bericht 6, (1968). [5] B. Ronnenberg, Ein selbstjustierendes 3-Komponenten-Laserdoppleranemometer nach dem Vergleichsstralverfahren, angewandt für Untersuchungen in einer stationären zylinder-symmnetrischen Drehströmung mit einem Rückstromgebiet, Max-Planck-Institut für Strömungsforschung, Göttingen, Bericht 20, (1977). [6] C.T. Hill, J.D. Huppler, R.B. Bird, Secondary flows in the disk-and-cylinder system, Chem. Eng. Sci. 21 (1966) 815–817. [7] M.P. Escudier, Observations of the flow produced in a cylindrical container by a rotating endwall, Exp. Fluids 2 (1984) 189–196. [8] H.J. Lugt, H.J. Haussling, Axisymmetric vortex breakdown in rotating fluid within a container, Trans. ASMEE: J. Appl. Mech. 49 (1982) 921–929. [9] H.J. Lugt, M. Abboud, Axisymmetric vortex breakdown with and without temperature effects in a container with a rotating lid, J. Fluid Mech. 179 (1987) 179–200. [10] J.N. Sorensen, High-order axisymmetric Navier-Stokes code: description and evaluation of boundary conditions, Int. J. Numer. Meth. Fluids 9 (1989) 1517–1537.

366

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

[11] J.M. Lopez, Axisymmetric vortex breakdown, Part 3. Onset of periodic flow and chaotic advection, J. Fluid Mech. 221 (1990) 533–552. [12] G.L. Brown, J.M. Lopez, Axisymmetric vortex breakdown. Part 2. Physical mechanisms, J. Fluid Mech. 221 (1990) 553– 576. [13] J.M. Lopez, A.D. Perry, Axisymmetric vortex breakdown. Part 3. Onset of periodic flow and chaotic advection, J. Fluid Mech. 234 (1992) 449–471. [14] A. Spohn, M. Mory, E.J. Hopfinger, Experiments on vortex breakdown in a confined flow generated by a rotating disc, J. Fluid Mech. 370 (1998) 73–99. [15] C.T. Hill, Nearly viscometric flow of viscoelastic fluids in the disk and cylinder system. II. Experimental, Trans. Soc. Rheol. 16(2) (1972) 213–245. [16] C. Day, J. Harris, I. Soria, D.V. Boger, M.C. Welsh, Behavior of an elastic fluid in cylindrical swirling flow, Exp. Therm. Fluid Sci. 12 (1996) 250–255. [17] M.P. Escudier, L.M. Cullen, Flow of a shear-tinning liquid in a cylindrical container with a rotating end wall, Exp. Therm. Fluid Sci. 12 (1996) 381–387. [18] D.V. Boger, A highly elastic constant-viscosity fluid, J. Non-Newtonian Fluid Mech. 3 (1977/78) 87–91. [19] J.R. Stokes. L.J.W. Graham, D.V. Boger, Observation of elastic effects in confined swirling flow of viscoelastic fluids, in: R.W. Bilger (Ed.), Proc. 12th Aust. Fluid Mech. Conf., University of Sydney, Sydney, Australia 1995, pp. 767–770. [20] J.R. Stokes, L.J.W. Graham, D.V. Boger, Vortex breakdown in confined swirling flow of a dilute flexible polymer solution, in: A. Ait-Kadi. J.M. Dealy, D.F. James, M.C. Williams (Eds.), Proc. XIIth Int. Congr. on Rheol. Laval University, Quebec City, Canada, 1996, pp. 359–360. [21] J.R. Stokes. L.J.W. Graham, D.V. Boger, Instabilities in the swirling flow of elastic liquids, in: Q.D. Nguyen, R.R. Huilgol, (Eds.), Proc. 8th National Conf. on Rheol., University of Adelaide, Adelaide, Australia, 1998, pp. 201–203. [22] J.M. Kramer, M.W. Johnson, Nearly viscometric flow in the disk and cylinder system. I. Theoretical, Trans. Soc. Rheol. 16(2) (1972) 197–212. [23] J.P. Nirschl, W.E. Stewart, Computation of viscoelastic flow in a cylindrical tank with a rotating lid, J. Non-Newtonian Fluid Mech. 16 (1984) 233–250. [24] F. Chiao, Hsueh-Chia Chang, Instability of a Criminale–Ericksen–Fibery fluid in a disk- and cylinder system, J. Non-Newtonian Fluid Mech. 36 (1990) 361–394. [25] O. Wunsch, G. Bohme, On torsionally driven viscoelastic flow in a cylindrical vessel, in: A. Ait-Kadi. J.M. Dealy, D.F. James, M.C. Williams (Eds.), Proc. XIIth Int. Congr. on Rheology, Laval University, Quebec City, Canada, 1996, pp. 359–360. [26] T. Moroi, M. Itoh, H. Toda, Numerical simulation of viscoelastic flow due to a rotating disc enclosed in a cylindrical casing, in: M.C. Thompson, K. Hourigan (Eds.), Proc. 13th Aust. Fluid Mech. Conf., Monash University, Melbourne, Australia 1998, pp. 313–316. [27] R.B. Bird, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids. Vol. 1. Fluid Mechanics, Wiley, New York, 1977. [28] S.-C. Xue, N. Phan-Thien, R.I. Tanner, Three-dimensional numerical simulations of viscoelastic flows through planar contractions, J. Non-Newtonian Fluid Mech. 74 (1998) 195–245. [29] S.-C. Xue, N. Phan-Thien, R.I. Tanner, Numerical investigations of Lagrangian unsteady extensional flows of viscoelastic fluids in 3-D rectangular ducts with sudden contractions, Rheol. Acta 37 (1998) 158–169. [30] N. Phan-Thien, R.I. Tanner, A new constitutive equation derived from network theory, J. Non-Newtonian Fluid Mech. 2 (1977) 353–365. [31] R.I. Tanner, in: Paper presented at the 6th Workshop on Numerical Methods, Denmark, 1989. [32] S.-C. Xue, Three-dimensional finite volume modelling and numerical simulations of viscoelastic flows, Ph.D. Thesis, The University of Sydney, Australia, 1997. [33] D. Rajagopalan, R.C. Armstrong, R.A. Brown, Finite element methods for calculation of steady viscoelastic flow using constitutive equations with a Newtonian viscosity, J. Non-Newtonian Fluid Mech. 36 (1990) 159–192. [34] J. Sun, N. Phan-Thien, R.I. Tanner, An adaptive viscoelastic stress splitting scheme and its application: AVSS/SI and AVSS/SUPG, J. Non-Newtonian Fluid Mech. 65 (1996) 75–91. [35] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York, 1980. [36] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Mech. Appl. Mech. Eng. 32 (1982) 199–259.

S.C. Xue et al. / J. Non-Newtonian Fluid Mech. 87 (1999) 337–367

367

[37] X.-L. Luo, R.I. Tanner, A decoupled finite element streamline-upwind scheme for viscoelastic flow problems, J. Non-Newtonian Fluid Mech. 31 (1989) 143–162. [38] X.-L. Luo, A control volume approach for integral viscoelastic model and its application to contraction flow of polymer melts, J. Non-Newtonian Fluid Mech. 64 (1996) 173–189. [39] C.M. Rhie, W.L. Chow, Numerical study of the turbulent flow past an airfoil with tailing edge separation, AIAA J. 21 (1983) 1525–1529. [40] J.P. Pulicani, J. Ouazzani, A Fourier–Chebyshev pseudospectral method for solving steady 3-D Navier–Stokes and heat equations in cylindrical cavities, Comput. Fluids 20(2) (1991) 93–109. [41] C. de Valh Davis, A note on a mesh for use with polar coordinates, Numer. Heat Transfer 2 (1979) 261–266. [42] L. Ma, D.B. Ingham, X. Wen, A finite volume method for fluid flow in polar cylindrical grids, Int. J. Numer. Meth. Fluids 28 (1998) 663–677. [43] H. Ozoe, K. Toh, A technique to circumvent a singularity at a radial center with application for a three-dimensional cylindrical system, Numer. Heat Transfer Part B 33 (1998) 355–365. [44] S.-C. Xue, N. Phan-Thien, R.I. Tanner, Numerical study of secondary flows of viscoelastic in straight pipes by an implicit finite volume method, J. Non-Newtonian Fluid Mech. 59 (1995) 191–213. [45] R.I. Issa, Solution of implicitly discretised fluid flow equations by operator-splitting, J. Comput. Phys. 62 (1985) 40–65. [46] J.A. Michelsen, Modeling of laminar incompressible rotating fluid flow, AFM 86-05. Ph.D. Dissertation, Department of Fluid Mechanics, Technical University of Denmark, 1986. [47] J.R. Stokes, private communication (1998). [48] H. Giesekus, A simple constitutive equation for polymer fluids based on the concept of deformation dependent tensorial mobility, J. Non-Newtonian Fluid Mech. 11 (1982) 69–109. [49] J.H. Ahlberg, E.N. Nilson, J.L. Walsh, The Theory of Splines and Their Applications, Academic Press, New York, 1967.