A locally-upwinded spectral technique (LUST) for viscoelastic flows

J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

A locally-upwinded spectral technique (LUST) for viscoelastic flows Robert G. Owens a,∗ , Cédric Chauvière a,1 , Timothy N. Philips b a

b

FSTI-ISE-LMF, Ecole Polytechnique Fédérale de Lausanne, CH 1015 Lausanne, Switzerland Department of Mathematics, University of Wales, Aberystwyth, SY23 3B2, Ceredigion, Wales, UK Received 25 July 2001; received in revised form 7 February 2002

Abstract In this paper, we develop an SUPG spectral element scheme suitable for computations of viscoelastic flows at high Deborah numbers. The novelty of the scheme lies in the derivation of the upwinding factors used in the perturbed test tensors of the weak form of the equations. These factors are related to the Deborah number and to the local mesh spacing and have been derived using the superconsistency ideas of Funaro [SIAM J. Numer. Anal. 30 (1993) 1664–1676; J. Sci. Comp. 12 (1997) 385–394; Comput. Math. Appl. 33 (1997) 95–103]. Results of applying the method to the benchmark problem of flow past a single confined cylinder demonstrate the superior stability and accuracy properties of the new scheme when compared with the SUPG spectral element method used previously by the authors [J. Non-Newtonian Fluid Mech. 95 (2000) 1–33; Comput. Meth. Appl. Mech. Engrg. 190 (2001) 3999–4018]. The new method gives excellent agreement with reference results in the literature. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Locally-upwinded spectral technique; Viscoelastic flows; SUPG

1. Introduction Spectral element methods were first introduced to the world some 18 years ago by Patera [7,35,41]. Together with their close cousin the p finite element method, spectral element methods are stable for linear elliptic and hyperbolic problems and exhibit exponential rates of convergence for smooth problems. We shall refer to these two methods together as high-order element methods in the discussion that follows. Given the greater computational challenge posed by the system of nonlinear partial differential equations of mixed mathematical type that govern certain viscoelastic flows, it comes as no surprise to note that high-order element methods have only been successfully employed for the solution of such systems of ∗ Corresponding author. Tel.: +41-21-6933589; fax: +41-21-6935307. E-mail address: [email protected] (R.G. Owens). 1 Present address: Division of Applied Mathematics, Brown University, 182 George Street, P.O. Box F, Providence, RI 02912, USA.

0377-0257/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 0 2 5 7 ( 0 2 ) 0 0 1 2 4 - 6

50

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

equations within the last decade. Since the rate of convergence of this class of methods is determined by the smoothness of the exact solution to the problem being solved, it has been natural to apply these methods to non-Newtonian flows in geometries devoid of singularities. Popular amongst these have been flows: • • • •

in periodically constricted tubes [29,48–50,53,54], through square arrays of cylinders [29,49,48], between eccentrically rotating cylinders (journal bearing problem) [25–27,31,32,42], past a single confined sphere or cylinder [11,12,16,39,40,56,57].

It should be added that some investigation of the performance of high-order element methods for viscoelastic flows with geometric singularities has been done and in particular for flow of fluids of the Oldroyd/Maxwell type through a sharp contraction [55,56] and for simplified one-dimensional models [49,50]. A review of papers featuring spectral collocation and pseudo-spectral methods for viscoelastic flows may be found in [12]. In the first two benchmark problems cited in the list above, the contributions of Khomami and co-workers feature prominently and are concerned with the use of hp Galerkin finite element methods for flows of Oldroyd-B and UCM fluids. In [29] and [49], hp Galerkin finite element methods were used for solving steady creeping flow through undulating tubes and through square arrays of cylinders at different levels of porosity. A comparison in terms of stability, accuracy and computational efficiency of the hp scheme with various lower-order finite element schemes (Galerkin, EVSS/SU, EVSS/SUPG) for flow of a UCM fluid revealed that only the hp, EVSS/SU and EVSS/SUPG methods were stable, and that moreover the hp scheme manifested exponential rates of convergence towards the exact solution, in contrast with the linear rates of the lower-order schemes. For steady, inertial flows of a UCM fluid in an undulating tube under change of type conditions Talwar and Khomami [50] found that an EVSS hp Galerkin scheme had an advantage over their previous hp scheme (no stress splitting) in that no Mach or elasticity number limitations were detected. Lower-order EVSS/SU methods for such problems, although convergent towards the exact solution, converged at an extremely slow rate, thus making their use for this class of problem prohibitively expensive. Viscoelastic flow (with and without inertia) through an undulating tube was the problem of choice of Van Kemenade and Deville in 1994 [53,54]. Their work represented the first application of Legendre spectral element methods to viscoelastic flows. In [53], the authors solved creeping flow of a UCM fluid through a corrugated tube. The limiting Weissenberg number encountered by the authors seemed to depend heavily upon the distribution of spectral elements. Nevertheless, excellent agreement for calculations of the flow resistance with the results of the pseudo-spectral and pseudo-spectral finite difference methods of Pilitsis and Beris [43,45] and of the 4 × 4 SUPG method of Marchal and Crochet [37] was demonstrated. No limiting Weissenberg number was reported for the finest spectral element meshes used, and the method was cheaper from the point of view of CPU time and degrees of freedom than the 4 × 4 SUPG method at comparable levels of accuracy. In a follow-up paper [54], the authors examined inertial flows and although good agreement with the flow resistance results from the PCFD pseudo-spectral method of Pilitsis and Beris [44] at a Reynolds number of 12 could be shown, the increase of the flow complexity due to inertia as the Weissenberg number increased led to a breakdown of the calculations at a lower Weissenberg number than had been achievable in [53]. The authors expressed their conviction that the numerical failure of their Galerkin spectral element method to reach high Weissenberg numbers in the case of flows with inertia did not trace its origin in the

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

51

order of the space discretization and that an investigation into some form of stabilization could prove useful. As further evidence of the successful application of Galerkin spectral element methods to non-Newtonian flows in smooth geometries, mention should be made of the series of papers from the Aberystwyth group on the journal bearing problem [25–27,31,32,42]. The great majority of these deal with generalized Newtonian fluids, albeit with viscosities that depend on pressure, temperature and shear-rate. However, in [42], and drawing inspiration from the viscosity law used earlier by Davies and Li [13], Phillips et al. used a Galerkin spectral element method and a White–Metzner model to investigate the performance of both statically and variably loaded journal bearings. In the simulations, the relaxation time/viscosity ratio was maintained constant for each run, but an assessment of the stability of the calculations is difficult since the ratios used were small (O(10−4 )) and no Weissenberg numbers were supplied. We turn from our optimism about the general applicability of Galerkin high-order element methods to viscoelastic flows at this point, however. Despite the successes showcased above, the fact remains that there are even geometrically smooth problems where the Galerkin formulation has proved to be a miserable failure when compared with low-order stabilized finite element or finite volume methods. Here we are thinking, for example, of viscoelastic flow past a single confined sphere or cylinder. Some early attempts to use Galerkin spectral element methods for these problems (see [40] and the Galerkin results presented in [11,12], for example) have been plagued by poor accuracy, low Deborah number limits and spurious numerical oscillations. Although these problems are defined in geometries that are smooth, huge stress wakes, thin elastic boundary layers, and a complex mixture of shear and extensional flow regions combine forces to ensure that the numerical solution of flow past a single confined sphere or cylinder is a formidable challenge, not be undertaken lightly by the computationally faint of heart. At the time of writing, convergence at Deborah numbers exceeding unity for high-order element simulations of the flow of a UCM or Oldroyd-B fluid past either a sphere in a tube or a cylinder in a channel at aspect ratios of 1/2, has not been possible with a classical mixed Galerkin method [11,12,16,56,57]. In 1995, Fan and Crochet [16] developed a high-order p finite element method and demonstrated p convergence of their numerical solution of UCM fluid flow in the 1:2 sphere-in-a-tube geometry up to a Deborah number of 2.1. The discretization of the constitutive equation was stabilized using the SUPG formulation. Denoting a finite-element length in the direction of flow by h, the perturbed test tensor S˜ in the SUPG formulation could be written: S˜ = S + huˆ · ∇S, (1) where S is the original test tensor and uˆ some non-dimensional scaling of the velocity field u. Three different possibilities of constructing uˆ were considered by Fan and Crochet: defining uˆ to be a unit vector in the direction of u (SUPG1), doing as did Lunsmann et al. [34] and dividing u by some characteristic velocity magnitude U (SUPG2) or dividing u in each element by an average taken over the nodes of the element (SUPG3). Best results were obtained with SUPG2 and SUPG3, and the importance of the upwinding factor huˆ vanishing on stationary walls for p convergence of the method was noted. In 1996–1997, Warichet and Legat [56,57] also undertook the solution of the 1:2 sphere-in-a-tube problem; this time with an adaptive high-order p finite element method. In [57], the SUPG2 method of Lunsmann et al. [34] was employed to stabilize the computations and this, in conjunction with an AVSS formulation, yielded converged solutions up to a Deborah number of 2.49 at much lower cost than with a uniform p AVSS/SUPG implementation of the method. It should be noted that in neither the work of Fan and Crochet nor of Warichet and Legat was any dependence of the upwinding factor huˆ upon the

52

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

polynomial order of representation of the dependent variables made. The length scale h was related only to the element diameter. Streamline upwinding was also found necessary for the sphere and cylinder benchmark problems by Chauvière and Owens [11,12] in their spectral element computations of flows of an Oldroyd-B fluid. The constitutive equation was solved on an element-by-element basis, the stress outflow information from an element being used as inflow data for the elements immediately downstream, and the SUPG2 formulation was adopted to further stabilize the scheme. The length scale h was chosen, in an ad hoc fashion to equal 1/N 2 in each spectral element, where N was the polynomial order in each spatial direction of the discrete representation of the components of velocity and stress. Error indicators developed in [11] showed that for the sphere problem the streamline upwinded spectral element method produced a solution that had a relative percentage error approximately a half of that of the Galerkin method when used on the same grid at De = 1.2. In the case of the cylinder problem, results in [12] demonstrated that the use of the element-by-element method and SUPG2 allowed limiting Deborah numbers more than double those possible with a standard Galerkin spectral element method to be obtained. This was partly due to the finer meshes that could now be afforded by solving the constitutive equation on an element-by-element basis but mainly due to the stabilization effect of SUPG2. It is the purpose of this paper to derive a new SUPG upwinding parameter huˆ and in particular to do so locally within each member of the partition of a spectral element formed by the quadrature grid. The derivation of the method draws inspiration from the operator-dependent grid ideas of Funaro [19–21] and from bubble stabilization methods [6,8,9]. We demonstrate why the proposed choice of upwinding factor appears reasonable and show that it is effective in stabilizing our calculations, whilst preserving good accuracy properties. The paper is organized as follows: in Section 2, we outline the problem to be solved by writing down the system of nonlinear coupled equations governing the flow of a rather general class of viscoelastic fluids whose constitutive equations are of differential form. Section 3 is dedicated to a discussion of appropriate choices of approximation spaces for the pressure, velocity, and elastic extra-stresses in order that the discrete Galerkin formulation of the equations should represent a well-posed problem. Then in Section 4, we describe a manner in which the Galerkin formulation of the constitutive equation may be modified by perturbing the test tensor. This is shown to lead to a consistent SUPG scheme. The upwinding factors in the test tensor perturbation are shown to be limited in magnitude by the dimensions of the mini-element formed from the partition of a spectral element by the Gauss–Lobatto grid, and an efficient algorithm for their computation is explained. The beneficial effects of locally upwinding the discrete equations is demonstrated in Section 5 by applying the method to the benchmark problem of a viscoelastic flow past a single confined cylinder. Finally, in Section 6 we draw some conclusions. 2. Problem description Let Ω ⊂ Rd be some bounded Lipschitz domain having boundary ∂Ω. We consider steady, inertialess, isothermal flow of a viscoelastic fluid having governing equations as follows: ∇ · u = 0,

(2)

−η2 ∇ u + ∇p = ∇ · τ + b,

(3)

2



τ + λτ + f (τ , γ˙ ) = η1 γ˙ (u),

(4)

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

53

where u denotes the fluid velocity, p the pressure, τ the elastic extra-stress tensor, b a body force, γ˙ (u) the rate-of-strain tensor and η ≡ η1 + η2 is the fluid viscosity. The tensor-valued function f of the elastic extra-stress tensor and the rate-of-strain tensor allows for possibly nonlinear departures from models of the Oldroyd–Maxwell type, such as the PTT model and Giesekus/Leonov models. The differential operator ∇ in (4) is the usual upper-convected derivative. 2.1. Weak formulation Let us denote the function spaces for the velocities, pressure and stresses by V , Q, and Σ, respectively. With V , Q, and Σ suitably defined (see, for example, [11] for a discussion on this matter in the case of models of the Oldroyd/Maxwell type) we may now define bilinear forms a : V ×V → R, b : Q×V → R, c : Σ × V → R, thus:
a(w, v) = η2 ∇w T : ∇vdx ∀w, v ∈ V , (5) Ω


b(q, v) =

Ω

q∇ · vdx

∀q ∈ Q, v ∈ V ,

(6)


c(T , v) =

Ω

T : ∇vdx

∀T ∈ Σ, v ∈ V .

(7)

We may also introduce a trilinear form d : V × Σ × Σ → R defined by:

 d(w, T , S) = (T + λT + f (T , γ˙ )) : S dx − η1 γ˙ (w) : S dx, Ω

Ω

∀w ∈ V , T , S ∈ Σ,

(8)

With a(·, ·), b(·, ·), c(·, ·) and d(·, ·, ·) defined as in (5)–(8) the Galerkin weak form of the governing equations in the absence of body forces is: find (u, p, τ ) ∈ V × Q × Σ such that: a(u, v) − b(p, v) + c(τ , v) = 0, b(q, u) = 0,

∀v ∈ V ,

(9)

∀q ∈ Q,

d(u, τ , Σ) = 0,

(10)

∀Σ ∈ Σ.

(11)

3. Discretization (Galerkin method) Let V N ⊂ V , QN ⊂ Q, and Σ N ⊂ Σ be conforming discrete subspaces of V , Q and Σ, respectively. Then a discrete formulation of (9)–(11) is: find (uN , pN , τ N ) ∈ V N × QN × Σ N such that: a(uN , v N ) − b(p N , v N ) + c(τ N , v N ) = 0, b(q N , uN ) = 0,

∀q N ∈ QN ,

d(uN , τ N , S N ) = 0,

∀S N ∈ Σ N .

∀v N ∈ V N ,

(12) (13) (14)

In this paper, we will use spectral element methods [7,35,41] for the discretization of the continuous problem (9)–(11). In [12], the nonlinear system of Eqs. (12)–(14) was solved in the case f ≡ 0 (Oldroyd-B

54

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

fluid) using Newton’s method with a finite difference approximation to the Jacobian and a preconditioned GMRES iterative method [46]. Because of the hyperbolic character of the constitutive Eq. (4), it proved possible in [12] to integrate it on an element-by-element basis after a careful ordering of the spectral elements in such a way that outflow stress data from elements further up the ordering was always available as inflow data for elements immediately downstream. Such an ordering always exists for flows without recirculation [30] and the reader is referred to [12] for further details. In this paper, we will use the same element-by-element approach and nonlinear solver as in [12]. An important issue in mixed spectral/finite element discretizations of the governing equations for Newtonian flows and viscoelastic flows is that of compatibility of the approximation spaces so as to ensure a well-posed discrete problem (12)–(14). On account of the intractability of the full set of equations, the available analyses (see, for example, [3,18,23]) are restricted to linearized versions of the constitutive equations for slow flows or the three fields Stokes problem obtained from (12)–(14) by setting λ = 0. From the work of Fortin and Pierre [18] and Baranger and Sandri [3] in the finite element context, and Gerritsma and Phillips [23] for spectral element methods, it may be shown that the usual LBB or inf-sup condition [5] on the spaces QN and V N applies. For spectral element methods one way of ensuring satisfaction of the LBB condition in either the primitive variable formulation of the Stokes problem [36] or the three fields Stokes problem [23] is to seek tensorized bases for the pressure consisting of polynomials in each spatial variable of degree 2 less in each element than those used for the velocity components. Lagrange interpolating polynomials based on, for example, the Gauss–Lobatto–Legendre (GLL) quadrature points [7] may be used for the velocity representation and the interior GLL points may be used in the pressure approximation. Let us partition the domain Ω into K (say) non-overlapping spectral elements {Ωk }K k=1 , and let Nk denote a d-tuple (Nk,1 , . . . , Nk,d ), for any 1 ≤ k ≤ K. We may then denote by PNk (Ωk ) the space of all polynomials on Ωk of degree ≤ Nk,j in the j th spatial direction, and further define: PN (Ω) = {φ : φ|Ωk ∈ PNk (Ωk )}.

(15)

The velocity and pressure approximation spaces may then be defined respectively as: V N = V ∩ (PN (Ω))d ,

(16)

QN = Q ∩ PN−2 (Ω),

(17)

where V N is a continuous approximation space and QN is discontinuous [11]. Of course, other choices of compatible discrete spaces exist [4]. In the mixed Galerkin formulation (12)–(14) of the three fields Stokes limit of the UCM equations (λ = 0, η2 = 0) a further compatibility condition applies between V N and the discrete stress space Σ N [18,23]. One way of ensuring compatibility is to choose the discrete rate-of-strain tensor to lie in Σ N and therefore to choose discontinuous stresses [23]. Some heuristic guidelines also exist in the literature for the choice of relative dimensions of V N and Σ N in the case of high-order mixed formulations for the full UCM equation [48,53,54]. Baranger and Sandri [3] showed, however, that no V N − Σ N compatibility condition is required for the Stokes limit of the Oldroyd-B equations (λ = 0, η2 = 0). For smooth problems, even in the case of a UCM fluid, there would appear to be no obvious advantage in taking the stresses to be discontinuous [24]. Noting this and exulting in the freedom that is ours, we choose the discrete stress space Σ N for our Oldroyd-B calculations described in Section 5 to be Σ N = {T ∈ Σ : Ti,j ∈ PN (Ω) ∩ C 0 (Ω)

∀i, j }.

(18)

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

55

That is, Σ N is a continuous space spanned by the same degree polynomials as used for the velocity components.

4. A LUST for viscoelastic flow calculations We begin our explanation of LUST by writing down a consistent SUPG scheme for the constitutive equation (4) with an upwinding factor h. Thus the weak form of the constitutive equation in LUST becomes d(u, τ , S + hu · ∇S) = 0,

∀S ∈ Σ.

(19)

The obvious question to ask now is how we define the scalar upwinding factor h. In fact, since we discretize Eq. (19) using a quadrature rule based on the GLL points in each element, all we need to know are the evaluations hij , where hij corresponds to a tensor S N ij all of whose components bar one are zero, the non-zero component being a Lagrange polynomial based on the GLL points and vanishing at all of them except for the (i, j )th. The upwinding factor is thus dependent upon the test tensor and not on spatial position. 4.1. Definition and properties of hij Let us denote by PN+1 (x) the non-zero polynomial of minimum degree (unique up to multiplication by a constant) that vanishes at all the GLL points in a given spectral element. From this point onwards and for ease of notation and clarity of explanation, we will work in Cartesian coordinates and consider only the two-dimensional parent element: Ωˆ = {(ξ, η) : −1 ≤ ξ, η ≤ 1}. The generalization to a physical element simply involves the introduction of a suitable mapping from the physical element to the parent element but the principles of the derivation of the upwinding factor remain the same. The generalization to three dimensions is entirely straightforward. It should also be noted that although the context in which the LUST is derived here is that of high-order element methods, analogous arguments should lead to an SUPG upwinding factor in the case of low-order finite element methods too. From the definition of the GLL quadrature points [7] PN+1 may be taken to be PN+1 (x) = PN+1 (ξ, η) = (1 − ξ 2 )LN (ξ )(1 − η2 )LN (η),

(20)

where LN is the degree N Legendre polynomial. Following the superconsistency ideas of Funaro [19–21] we than define hij so that N N d(uN , τ N + αΠ N+1 , S N ij + hij u · ∇S ij ) = 0,

(21)

where α is an arbitrary constant and Π N+1 is the second-order tensor each of whose components is the polynomial PN+1 . Interpreting (21) we see that our choice of upwinding factor hij allows the discrete stress trial space Σ N to be augmented by the addition of a linear space spanned by Π N+1 . For the sake of simplicity, we will work from this point onwards with fluids of the Oldroyd/Maxwell type. In this case,

56

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

Fig. 1. Parent element Ωˆ = [−1, 1]2 partitioned by the GLL grid. Ωˆ ij is a typical mini-element.

the differential operator in (4) is linear in the stress tensor and using integration by parts the discretized equations for LUST become N N d(uN , τ N , S N ij + hij u · ∇S ij )GLL = 0,

(22)

with hij chosen so that [L(Π N+1 ) − hij uN · ∇L(Π N+1 )](x ij ) = 0,

(23)

where d(·, ·, ·)GLL denotes the GLL quadrature evaluation of d(·, ·, ·) and L is differential operator in (4). We note that (23) is the linear part of the Taylor series expansion of L(Π N+1 )(zij ) = 0,

(24)

where zij = x ij −hij uN (x ij ). Therefore, assuming hij uN (x ij ) to be sufficiently small and neglecting velocity gradients by considering that the components of velocity are constants in each of the “mini-elements” ˆ ij ⊂ ˆ formed from the GLL grid (see Fig. 1), (24) may be replaced, for each component of the tensor, by the approximate scalar equation PN+1 (zij ) + λ(uN · ∇)PN+1 (zij ) = 0.

(25)

Since, for λ = 0 the solutions, by definition of PN+1 , are just the points {x ij } of the GLL grid, it is natural to associate each root zij of (25) with the point x ij to which it would collapse as λ → 0. In Fig. 2, we consider a mini-element having one of its four vertices the GLL point x ij and suppose that at that point the velocity vector uN is directed towards the exterior of the mini-element. We term such a vertex an outflow vertex for this mini-element, for obvious reasons. We further suppose, with no loss of generality, that PN+1 ≥ 0 in this mini-element. At x ij it is clear that PN+1 − λ(uN · ∇)PN+1 ≤ 0, since PN+1 vanishes at x ij and the streamline derivative must be non-positive. Now suppose that λuN (x ij ) = 0. Then, backtracking along the streamline Γ (x ij ) passing through x ij (see Fig. 2) we see that there must

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

57

Fig. 2. Diagram showing the positioning of the zero zij of PN +1 − λu · ∇PN +1 in a mini-element having an outflow vertex x ij .

exist a point x ∗ij (say) where PN+1 attains a maximum along the streamline and therefore (uN · ∇)PN+1 (x ∗ij ) = 0. Hence, at x ∗ij PN+1 + λ(u · ∇)PN+1 ≥ 0,

(26)

and by the intermediate value theorem there must exist a zero zij lying on the streamline between x ij and x ∗ij . If Γ (x ij ) runs along the edge of the mini-element or if λ = 0, then zij = x ij . Otherwise, zij lies in the interior of the mini-element for which x ij is an outflow vertex. We see that the evaluation of (22) using a GLL quadrature rule requires us only to find the shifts {hij } corresponding to the test tensor S N ij . We describe in the next subsection how to do this. 4.2. Computation of the shift factors hij In order to compute the shift factor hij corresponding to a test tensor S N ij , we may expand PN+1 (zij ) as a Taylor series about x ij and truncate at some order in hij . Doing the same for λ(uN · ∇)PN+1 (zij ) and adding the expressions leads to PN+1 (zij ) + λ(u · ∇)PN+1 (zij ) N

h2ij

N

h3ij

= PN+1 (x ij ) − hij (u · ∇)PN+1 (x ij ) + (u · ∇) PN+1 (x ij ) − (uN · ∇)3 PN+1 (x ij ) + · · · 2 3!  2 hij + λ (uN · ∇)PN+1 (x ij ) − hij (uN · ∇)2 PN+1 (x ij ) + (uN · ∇)3 PN+1 (x ij ) 2  h3ij N − (u · ∇)4 PN+1 (x ij ) + · · · = 0. (27) 3 2

It is clear in (27) that PN+1 (x ij ) = 0. For points x ij not lying on the parent element boundary it may be shown, using the properties of Legendre polynomials and of the polynomial PN+1 , that the first

58

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

and third-order streamline derivatives appearing in (27) also vanish when evaluated at x ij . Hence, if we truncate the two Taylor series appearing in (27) at O(h3ij ), exploit the properties of derivatives of Legendre polynomials in order to obtain a common factor LN (ξi )LN (ηj ) in the equation for hij and finally divide (27) throughout by hij we get, after some rearrangement, the quadratic equation   u22,ij u21,ij 2λN(N + 1) + h2ij + hij − 2λ = 0. (28) 3 (1 − ξi2 ) (1 − ηj2 ) One of the roots of this equation is positive and this is what is chosen for hij . Of course, the decision to truncate the Taylor series in (27) at O(h3ij ) is somewhat empirical. It would be interesting to see what happened if higher-order terms were retained. The only complication then would be that hij would be more difficult to solve for. Truncating (27) at lower orders does not lead to meaningful results. A similar analysis to that conducted above may be performed for the interior edge GLL points and the four corner GLL points. We will now consider the application of LUST to a classical benchmark problem.

5. Numerical example: viscoelastic flow past a cylinder in a channel We consider the flow of an Oldroyd-B fluid past a cylinder placed symmetrically in a channel (see Fig. 3). Despite the smoothness of the geometry, the problem is acknowledged to be exceedingly difficult to solve due to the development of both severe stress boundary layers on the cylinder surface and an elastic wake as the Deborah number increases. This problem has proved to be numerically even more difficult than the related sphere problem because for the same aspect ratio the planar flow past a cylinder in a channel undergoes stronger contraction and expansion than axisymmetric flow past a sphere in a tube. A wide range of numerical results exist in the literature for this problem and this has allowed us to assess the performance of our scheme. Liu et al. [33] used numerical simulation to investigate flow of polymer solutions around a periodic linear array of cylinders where the fluid was modelled by the Giesekus, FENE-P and Chilcott–Rallison (CR) models. The Oldroyd-B model was considered as a limiting case of the CR model as the maximum dumbbell length L was allowed to go to infinity. The authors used two different finite element formulations for solving the problem: the EVSS-G and the DEVSS-G formulations, the ‘G’ indicating the fact that the velocity gradient ∇u was approximated by a (continuous) variable G separate from the other fields. Liu et al. [33] succeeded in getting converged solutions up to a Deborah number of 1.05 for the problem of a single cylinder and aspect ratio 1/2. The meshes used required up to 67,495 degrees of freedom.

Fig. 3. Cylinder radius R placed symmetrically in a two-dimensional channel of half width H .

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

59

A new stabilized formulation featuring a stabilization term in the momentum equations where this consisted of a square residual of the continuity equation, was introduced by Fan et al. [17] in 1999 and presented a challenge to the widely accepted belief that the momentum equation should be made explicitly elliptic. The constitutive equation (UCM and Oldroyd-B) was handled with a consistent SUPG method, and hp finite elements were used for the discretization of the resulting set of weak equations. The resulting Galerkin least squares formulation (called MIX1 by the authors) was then tested on three benchmark problems: flow of a UCM fluid between eccentric cylinders, flow of a UCM fluid around a sphere in a tube and finally, flow of a UCM and Oldroyd-B fluid around a cylinder in a channel with aspect ratio 1/2. For the cylinder problem and using a UCM fluid comparisons were made with two other finite element formulations: MIX0, in which a Galerkin method was used for the momentum and continuity equations and the SUPG technique for the constitutive equation, and the DEVSS method. Fan et al. [17] defined a Deborah number in identical fashion to the Weissenberg number of Liu et al. [33]. For the UCM fluid, limiting Deborah numbers of 0.75 and 0.8 were encountered for MIX1 and DEVSS, respectively; the major difficulty in obtaining convergence being identified as the accurate determination of the axial normal stress profile in the wake of the cylinder. Convergence could not be obtained for any Deborah number using MIX0. For flow of an Oldroyd-B fluid (with solvent to total viscosity ratio 0.59) through the same geometry the hp MIX0 method performed as well as the MIX1 and DEVSS formulations but convergence was lost after a Deborah number of 1.05. Sun et al. [47] developed a numerical method with several different components: a combination of the AVSS and DEVSS methods for preserving the elliptic character of the momentum and continuity equation pair, a continuous approximation G for the velocity gradient field, and a discontinuous Galerkin method. The resulting DAVSS-G/DG method was used for solving the system of equations for the Oldroyd-B and Giesekus models and applied to the cylinder problem for aspect ratios of 1/2 and 1/8. Drag calculation results were compared for the Oldroyd-B fluid with those from the DEVSS-G method of Liu et al. [33] for the larger of the two ratios. The DAVSS-G/DG Oldroyd-B calculations converged for the narrow channel up to De = 1.85 and significantly, no axial normal stress profiles in the wake of the cylinder were presented by the authors. Hence no assessment of whether the authors achieved convergence in the critical region in the wake of the cylinder can be made. Drag calculation results agreed well with the DEVSS-G calculations of Liu et al. [33] up to a Deborah number of 1.0. Two further solutions for the benchmark problem of viscoelastic flow past a confined cylinder have appeared very recently in the literature [1,10]. In [1], Alves et al. used the high resolution finite volume method of Oliveira et al. [38]. Two high resolution schemes were implemented and assessed for the treatment of the convective terms in the Oldroyd/Maxwell constitutive equations: the MINMOD method of Harten [28] and the SMART method of Gaskell and Lau [22]. The discretized sets of equations for each variable were solved sequentially in a decoupled manner and conservation of mass and momentum by the pressure and velocity fields was ensured by use of the SIMPLEC algorithm [52]. The calculated drag coefficients for both flows of a UCM and Oldroyd-B fluid were in particularly good agreement with the results of Fan et al. [17] and all the more so after the use of a Richardson’s extrapolation. In [10], Caola et al. used a DEVSS-G/DG finite element discretization for solving the flow of an Oldroyd-B fluid past a single confined cylinder. The method involved an operator-splitting time integration method for the decoupling of the pressure and velocity field from the integration of the constitutive equation and the resulting generalized Stokes problem was solved using a BiCGStab iterative method [51]. Up to 751,100 degrees of freedom were used and a limiting Deborah number of 1.1 was encountered for the 2:1 cylinder problem.

60

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

Table 1 Total number of degrees of freedom and number of quadrature points for different values of the polynomial degree N N

Total dof

Number of nodes

12 13 14 15 16 17 18 19 20

17575 20590 23845 27340 31075 35050 39265 43720 48415

3031 3542 4093 4684 5315 5986 6697 7448 8239

Given the wide range of numerical results in the literature for the cylinder-in-a-channel problem when the aspect ratio Λ = R/H is taken equal to 0.5 and the solvent to total viscosity ratio η2 /η is taken equal to 0.59, we chose the same parameters for our numerical investigations. The computational domain extends a distance 20R upstream and downstream of the cylinder to ensure that the assumption of fully developed flow conditions at entry and exit is valid. We imposed no-slip conditions on the cylinder surface (x 2 + y 2 = R 2 ) and on the channel wall (y = H ). In order to save on computational cost we have assumed that the flow has a plane of symmetry along y = 0 so that only half of the domain needs to be considered. The flow domain is divided into 20 conforming spectral elements and polynomial degrees of up to Ni = (20, 20) (1 ≤ i ≤ 20) are used. We denote, for simplicity, the polynomial order of approximation in either direction of the parent element by N. Fig. 4 shows three spectral element grids for N = 16, 18, and 20 from top to bottom, respectively. The bold lines indicate the boundaries of the spectral elements and the intersection of the thin lines are the positions of the quadrature points x ij . The total number of degrees of freedom and the corresponding number of quadrature points is given in Table 1 for the whole range of polynomial degrees to be used in this paper (12 ≤ N ≤ 20). 5.1. Numerical results For all the results presented in this section the elastic extra-stress is computed using the element-byelement technique described in detail in [12]. In conjunction with the element-by-element technique, three different approaches will be tested: the first one is a standard mixed Galerkin method, the second is the SUPG method with a constant upwinding factor equal to 1/N 2 as used in [11,12], and the last uses the LUST described in the previous section. Starting at a Deborah number De = 0.5, the calculations are performed at Deborah numbers that are successively increased in steps of 0.1 for De ≤ 1.1 and of 0.05 for De > 1.1, until convergence is lost. The computations are performed for meshes having corresponding total numbers of degrees of freedom as reported in Table 1. Fig. 5 gives the maximum Deborah number achievable for the three methods for polynomial degrees ranging from N = 12 to 20. As should be expected, the Galerkin method lacks robustness since the maximum Deborah number achievable is limited to De = 0.8. No clear trend in the limiting Deborah number as a function of N is discernible and the maximum Deborah number achievable is strongly dependent on whether N is odd or even. The SUPG method of [11,12] pushes the maximum Deborah

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

61

Fig. 4. Twenty spectral elements for three levels of discretization (N = 16, 18, and 20).

number up to De = 1.1 and with LUST we go up to De = 1.2. Apart from a couple of small “accidents” at N = 18,19, the overall trend in the LUST calculations is that the limiting Deborah number increases with mesh refinement. The robustness of the LUST is therefore demonstrated. As we will see in the paragraphs that follow, this gain of robustness does not appear to have been achieved at the expense of accuracy. The most common quantity that may be found in the literature for the comparison of numerical results is the drag on the cylinder, F which may be expressed as follows:      
π  ∂uy ∂ux ∂ux F = + −p + 2η2 (29) + τxx cos θ + η2 + τxy sin θ Rdθ. ∂x ∂x ∂y 0

62

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

Fig. 5. Maximum Deborah number achievable for three methods and different levels of discretization (N = 12–20). Oldroyd-B fluid.

As has been pointed out by several authors [2,11,12,15–17], drags on cylinders or spheres are not, however, a good indicator of the global quality of the solution: the computed drag may not be sensitive to numerical oscillations that may appear at the interface between spectral elements, or to inaccuracies in the body wake, for example. However, an inaccurate computation of the drag is almost a guarantee that something is badly wrong in the critical part of the flow: the wake. Figs. 6–8 represent the value of the drag for different Deborah numbers and three different meshes N = 16, 18, and 20. Fig. 6 shows the Galerkin results, Fig. 7 the results of the SUPG method of [11,12], and Fig. 8 the LUST results. For all methods, convergence of the drag with mesh refinement is excellent at De = 0.5 and still good up to De = 0.7. However, convergence with mesh refinement beyond that point is not evident except for the LUST. In Fig. 9, where the contour of the xx component of the extra-stress is represented for the three methods for N = 16 and De = 0.5, we can see that the solutions are essentially wiggle-free and, at this Deborah number, the three methods are more or less equivalent. The scattering of the computed drag values at high Deborah numbers in Figs. 6 and 7 is a strong indication that poor accuracy should be anticipated in the flow domain. This suspicion is confirmed in Fig. 10, where for the same mesh as the one used in Fig. 9, we have represented the contours of τxx , but now at De = 0.8. Very large oscillations are visible in the wake of the cylinder at the interface between the spectral elements for (a) the Galerkin method and (b) the SUPG method of [11,12], whereas the LUST gives a wiggle-free contour plot over the entire domain. Even at a higher Deborah number, the contours remain smooth as can be seen in Fig. 11 which represents τxx at De = 1.2 for the finest mesh N = 20. In Table 2, we compare and contrast, for two levels of discretization (N = 12 and 16), the number of (inner) GMRES iterations and (outer) Newton steps required by our Galerkin method, SUPG method and LUST to converge for Deborah numbers in the range [0.7, 1.2]. The convergence criterion applied

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

63

Fig. 6. Drag against Deborah number for three levels of discretization (N = 16, 18, and 20). Galerkin method. Oldroyd-B fluid.

Fig. 7. Drag against Deborah number for three levels of discretization (N = 16, 18, and 20). SUPG method (h = 1/N 2 ). Oldroyd-B fluid.

64

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

Fig. 8. Drag against Deborah number for three levels of discretization (N = 16, 18, and 20). LUST. Oldroyd-B fluid. Table 2 Total number of GMRES iterations and, in brackets, number of Newton steps for two levels of discretization (N = 12 and 16) for the three methods investigated in this paper De

0.7 0.8 0.9 1.0 1.1 1.2

N = 12

N = 16

Galerkin

SUPG

LUST

Galerkin

SUPG

LUST

594 (5) 831 (6) – – – –

409 (4) 563 (5) – – – –

347 (4) 402 (4) 549 (5) – – –

674 (5) 904 (6) – – – –

552 (5) 693 (5) 977 (6) 1059 (6) 1609 (8) –

412 (4) 575 (5) 669 (5) 915 (6) 1023 (6) 1145 (6)

here is that the difference in value of all the components τij (1 ≤ i, j ≤ 2) of the polymeric extra-stress tensor at the nth and (n + 1)th Newton step should satisfy |τijn+1 (x k/ ) − τijn (x k/ )| < 5 × 10−5 , at all the GLL points x k/ throughout the flow domain Ω. The inner iterations are stopped once the Euclidean norm of the residual drops below 1 × 10−9 . As is to be expected, for all three methods the total number of GMRES iterations and Newton steps required for convergence increases with the Deborah number. However, it is also clear from the table that convergence is most rapid in the case of the LUST.

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

65

Fig. 9. Contours of τxx computed with a uniform mesh N = 16 at De = 0.5. Oldroyd-B fluid.

Table 3 gives the values of the drag computed with LUST and the three finest meshes N = 18, 19, and 20 as well as the results of Caola et al. [10] and Fan et al. (MIX0) [17]. The results presented in this table seem to be in very close agreement. Let us now take a look at the comparison between our results and those of others in the literature. This is done in Fig. 12. There, we have additionally included the results as supplied in the original papers of Dou and Phan-Thien (DAVSS-ω) [14], Sun et al. [47] and Alves et al. (SMART) [1], as well as the results (as supplied by Sun et al. [47]) of Liu et al. (DEVSS-G) [33].

66

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

Fig. 10. Contours of τxx computed with a uniform mesh N = 16 at De = 0.8. Oldroyd-B fluid.

The scale of the ordinate axis is chosen to emphasis any small differences between the results. The following comments may be made: 1. The only results that match well up to De = 0.7 are those of Alves et al., Fan et al., Caola et al. and ours. 2. For 0.7 ≤ De ≤ 1.0, two trends are visible: the results of Alves et al. and Fan et al., which match extremely well and slightly below, the results of Caola et al. with ours.

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

67

Fig. 11. Contours of τxx computed with a uniform mesh N = 20 at De = 1.2. LUST. Oldroyd-B fluid.

Table 3 Drag factor computed on uniform meshes (N = 18–20) with LUST. Also included are results of Caola et al. [10] and results of Fan et al. [17], Oldroyd-B fluid De

0.0 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2

N 18

19

20

132.357 118.827 117.766 117.268 117.214 117.519 118.123 118.976 –

132.357 118.828 117.772 117.284 117.250 117.586 118.247 119.197 –

132.357 118.827 117.775 117.291 117.237 117.503 118.030 118.786 119.764

Caola et al. (dof = 751,110)

Fan et al.

132.384 118.763 – – – – 117.783 118.031 –

132.36 118.83 117.78 117.32 117.36 117.80 118.49 – –

Our maximum Deborah number, De = 1.2 competes well with the results of Liu et al. [33] and Fan et al. [17] who both encountered limiting Deborah numbers of 1.05, of Alves et al. [1] who encountered a limiting Deborah numbers of 1.0 and Caola et al. [10] who encountered a limiting Deborah number of 1.1. In terms of the number of degrees of freedom, the LUST is also very competitive. Our finest mesh (N = 20) uses only 48,415 degrees of freedom and gives results close to those of Caola et al. who used 751,110 degrees of freedom on their finest mesh. Unfortunately, when it comes to considering the xx component of the elastic extra-stress in the wake of the cylinder, no convergent trend can be observed for Deborah numbers greater than or equal to 0.8. This is illustrated in Figs. 13 and 14 for three levels of discretization (N = 16, 18, and 20) with LUST. At a Deborah number of 0.7 (Fig. 13), the convergence is clear both on the cylinder surface and in the wake of the cylinder. At De = 1.1, (Fig. 14), convergence on the cylinder surface is still good and explains why there is still a convergent trend for the drag coefficient. Huge differences can be seen in the wake of the cylinder at this Deborah number, however. Similarly, inspecting the axial normal stress in the wake of the cylinder, Fan et al. [17] concluded (and we concur) that numerical solutions for the Oldroyd-B fluid at Deborah numbers higher than about 0.8 are probably numerical artifacts.

68

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

Fig. 12. Comparison of our drag (N = 20 LUST) with [1,10,14,17,33,47] Oldroyd-B fluid.

Fig. 13. Profiles of τxx on the cylinder surface and in the wake of the cylinder for different meshes at De = 0.7. Oldroyd-B fluid.

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

69

Fig. 14. Profiles of τxx on the cylinder surface and in the wake of the cylinder for different meshes at De = 1.1. Oldroyd-B fluid.

6. Conclusions In this paper, we have developed a new locally-upwinded spectral technique (LUST) for the stable integration of the governing differential equations for certain classes of viscoelastic flows. The important ingredient in the method is the manner in which the SUPG upwinding factors are defined. The extension of the ideas here to three dimensions and lower-order finite elements is straightforward. The LUST has been seen to be extremely competitive with other reference formulations when applied to the benchmark problem of the flow of an Oldroyd-B fluid past a single confined cylinder. References [1] M.A. Alves, F.T. Pinho, P.J. Oliveira, The flow of viscoelastic fluids past a cylinder: finite-volume high-resolution methods, J. Non-Newtonian Fluid Mech. 97 (2001) 207–232. [2] F.P.T. Baaijens, S.H.A. Selen, H.P.W. Baaijens, G.W.M. Peters, H.E.H. Meijer, Viscoelastic flow past a confined cylinder of a low density polyethylene melt, J. Non-Newtonian Fluid Mech. 68 (1997) 173–203. [3] J. Baranger, D. Sandri, A formulation of Stokes’s problem and the linear elasticity equations suggested by the Oldroyd model for viscoelastic flow, Math. Model Num. Anal. 26 (1992) 331–345. [4] C. Bernardi, Y. Maday, Uniform inf-sup conditions for the spectral discretization of the Stokes problem, Math. Models Methods Appl. Sci. 9 (1999) 395–414. [5] F. Brezzi, On the existence, uniqueness and approximation of saddle point problems arising from Lagrangian multipliers, RAIRO Anal. Numér. 8 (1974) 129–151. [6] C. Canuto, Stabilization of spectral methods by finite element bubble functions, in: C. Bernardi, Y. Maday (Eds.), Proceedings of ICOSAHOM.92, Montpellier 1992, North-Holland, Amsterdam, 1994, Comput. Methods Appl. Mech. Engrg. 116 (1994) 13–26.

70

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

[7] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, Berlin, Heidelberg, 1988. [8] C. Canuto, A. Russo, V. Van Kemenade, Stabilized spectral methods for the Navier–Stokes equations: residual-free bubbles and preconditioning, Comput. Methods Appl. Mech. Engrg. 166 (1998) 65–83. [9] C. Canuto, V. Van Kemenade, Bubble-stabilized spectral methods for the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 135 (1996) 35–61. [10] A.E. Caola, Y.L. Joo, R.C. Armstrong, R.A. Brown, Highly parallel time integration of viscoelastic flows, J. Non-Newtonian Fluid Mech. 100 (2001) 191–216. [11] C. Chauvière, R.G. Owens, How accurate is your solution? Error indicators for viscoelastic flow calculations, J. Non-Newtonian Fluid Mech. 95 (2000) 1–33. [12] C. Chauvière, R.G. Owens, A new spectral element method for the reliable computation of viscoelastic flow, Comput. Methods Appl. Mech. Engrg. 190 (2001) 3999–4018. [13] A.R. Davies, X.K. Li, Numerical modelling of pressure and temperature effects in viscoelastic flow between eccentrically rotating cylinders, J. Non-Newtonian Fluid Mech. 54 (1994) 331–350. [14] H.-S. Dou, N. Phan-Thien, The flow of an Oldroyd-B fluid past a cylinder in a channel: adaptive viscosity vorticity (DAVSS-ω) formulation, J. Non-Newtonian Fluid Mech. 87 (1999) 47–73. [15] Y.R. Fan, A comparative study of the discontinuous Galerkin and continuous SUPG finite element methods for computation of viscoelastic flows, Comput. Methods Appl. Mech. Engrg. 141 (1997) 47–65. [16] Y. Fan, M.J. Crochet, High order finite element methods for steady viscoelastic flows, J. Non-Newtonian Fluid Mech. 57 (1995) 283–311. [17] Y.R. Fan, R.I. Tanner, N. Phan-Thien, Galerkin/least-squares finite-element methods for steady viscoelastic flows, J. Non-Newtonian Fluid Mech. 84 (1999) 233–256. [18] M. Fortin, R. Pierre, On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows, Comput. Methods Appl. Mech. Engrg. 73 (1989) 341–350. [19] D. Funaro, A new scheme for the approximation of advection-diffusion equations by collocation, SIAM J. Num. Anal. 30 (1993) 1664–1676. [20] D. Funaro, Improving the performance of implicit schemes for hyperbolic equations, J. Sci. Comput. 12 (1997) 385–394. [21] D. Funaro, Some remarks about the collocation method on a modified Legendre grid, Comput. Math. Appl. 33 (1997) 95–103. [22] P.H. Gaskell, A.K.C. Lau, Curvature-compensated convective-transport—SMART, a new boundedness-preserving transport algorithm, Int. J. Num. Methods Fluids 8 (1988) 617–641. [23] M.I. Gerritsma, T.N. Phillips, Compatible spectral approximations for the velocity–pressure–stress formulation of the Stokes problem, SIAM J. Sci. Comput. 20 (1999) 1530–1550. [24] M.I. Gerritsma, T.N. Phillips, On the use of characteristic variables in viscoelastic flow problems, IMA J. Appl. Math. 66 (2001) 127–147. [25] D.Rh. Gwynllyw, T.N. Phillips, Preconditioned iterative methods for unsteady non-Newtonian flow between eccentrically rotating cylinders, SIAM J. Sci. Comput. 17 (1996) 1369–1394. [26] D.Rh. Gwynllyw, A.R. Davies, T.N. Phillips, A moving spectral element approach to the dynamically loaded journal bearing problem, J. Comput. Phys. 123 (1996) 476–494. [27] D.Rh. Gwynllyw, A.R. Davies, T.N. Phillips, On the effects of a piezoviscous lubricant on the dynamics of a journal bearing, J. Rheol. 40 (1996) 1239–1266. [28] A. Harten, High-resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49 (1983) 357–393. [29] B. Khomami, K.K. Talwar, H.K. Ganpule, A comparative study of higher- and lower-order finite element techniques for computation of viscoelastic flows, J. Rheol. 38 (1994) 255–289. [30] P. Lesaint, P.A. Raviart, On a finite element method for solving the neutron transport equation, in: C. de Boor (Ed.), Mathematical Aspects of Finite Elements in Partial Differential Equations, Academic Press, New York, 1974, pp. 89–123. [31] X.K. Li, D.Rh. Gwynllyw, A.R. Davies, T.N. Phillips, Three-dimensional effects in dynamically loaded journal bearings, Int. J. Num. Methods Fluids 29 (1999) 311–341. [32] X.K. Li, D.Rh. Gwynllyw, A.R. Davies, T.N. Phillips, On the influence of lubricant properties on the dynamics of two-dimensional journal bearings, J. Non-Newtonian Fluid Mech. 93 (2000) 29–59. [33] A.W. Liu, D.E. Bornside, R.C. Armstrong, R.A. Brown, Viscoelastic flow of polymer solutions around a periodic, linear array of cylinders: comparisons of predictions for microstructure and flow fields, J. Non-Newtonian Fluid Mech. 77 (1998) 153–190.

R.G. Owens et al. / J. Non-Newtonian Fluid Mech. 108 (2002) 49–71

71

[34] W.J. Lunsmann, L. Genieser, R.C. Armstrong, R.A. Brown, Finite element analysis of steady viscoelastic flow around a sphere in a tube: calculations with constant viscosity models, J. Non-Newtonian Fluid Mech. 48 (1993) 63–99. [35] Y. Maday, A.T. Patera, Spectral element methods for the incompressible Navier–Stokes equations, in: A.K. Noor, J.T. Oden (Eds.), State of the Art Surveys in Computational Mechanics, American Society of Mechanical Engineers, New York, 1989, pp. 71–143. [36] Y. Maday, A.T. Patera, E.M. Rønquist, The PN − PN −2 Method for the Approximation of the Stokes Problem, Publications du Laboratoire d’Analyse Numérique R92025, Université Pierre et Marie Curie, 1992. [37] J.M. Marchal, M.J. Crochet, A new mixed finite element method for calculating viscoelastic flow, J. Non-Newtonian Fluid Mech. 26 (1987) 77–114. [38] P.J. Oliveira, F.T. Pinho, G.A. Pinto, Numerical simulation of non-linear elastic flows with a general collocated finite-volume method, J. Non-Newtonian Fluid Mech. 79 (1998) 1–43. [39] R.G. Owens, A posteriori error estimates for spectral element solutions to viscoelastic flow problems, Comput. Methods Appl. Mech. Engrg. 164 (1998) 375–395. [40] R.G. Owens, T.N. Phillips, Steady viscoelastic flow past a sphere using spectral elements, Int. J. Num. Methods Engrg. 39 (1996) 1517–1534. [41] A.T. Patera, A spectral element method for fluid dynamics: laminar flow in a channel expansion, J. Comput. Phys. 54 (1984) 468–488. [42] T.N. Phillips, R.E. Need, A.R. Davies, B.P. Williamson, L.E. Scales, The Effect of Viscoelasticity on the Performance of Dynamically Loaded Journal Bearings, Technical Report 982639, SAE Technical Paper Series, 1998. [43] S. Pilitsis, A.N. Beris, Calculations of steady-state viscoelastic flow in an undulating tube, J. Non-Newtonian Fluid Mech. 31 (1989) 231–287. [44] S. Pilitsis, A.N. Beris, Viscoelastic flow in an undulating tube. II. Effects of high elasticity, large amplitude of undulation and inertia, J. Non-Newtonian Fluid Mech. 39 (1991) 375–405. [45] S. Pilitsis, A.N. Beris, Pseudospectral calculations of viscoelastic flow in a periodically constricted tube, Comput. Methods Appl. Mech. Engrg. 98 (1992) 307–328. [46] Y. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869. [47] J. Sun, M.D. Smith, R.C. Armstrong, R.A. Brown, Finite element method for viscoelastic flows based on the discrete adaptive viscoelastic stress splitting and the discontinuous Galerkin method DAVSS-G/DG, J. Non-Newtonian Fluid Mech. 86 (1999) 281–307. [48] K.K. Talwar, H.K. Ganpule, B. Khomami, A note on the selection of spaces in computation of viscoelastic flows using the hp-finite element method, J. Non-Newtonian Fluid Mech. 52 (1994) 293–307. [49] K.K. Talwar, B. Khomami, Application of higher order finite element methods to viscoelastic flow problems in porous media, J. Rheol. 36 (1992) 1377–1416. [50] K.K. Talwar, B. Khomami, Higher order finite element techniques for viscoelastic flow problems with change of type and singularities, J. Non-Newtonian Fluid Mech. 59 (1995) 49–72. [51] H.A. van der Vorst, Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 13 (2) (1992) 631–644. [52] J.P. Van Doormal, G.D. Raithby, Enhancements of the SIMPLE method for predicting incompressible fluid flows, Num. Heat Transfer 7 (1984) 147–163. [53] V. Van Kemenade, M.O. Deville, Application of spectral elements to viscoelastic creeping flows, J. Non-Newtonian Fluid Mech. 51 (1994) 277–308. [54] V. Van Kemenade, M.O. Deville, Spectral elements for viscoelastic flows with change of type, J. Rheol. 38 (1994) 291–307. [55] R.G.M. van Os, M.I. Gerritsma, A variable order spectral element scheme applied to the VPTS formulation of the Upper Convected Maxwell model, J. Non-Newtonian Fluid Mech. 108/1–3 (2002) 73–97. [56] V. Warichet, V. Legat, Adaptive hp-finite element viscoelastic flow calculations, Comput. Methods Appl. Mech. Engrg. 136 (1996) 93–110. [57] V. Warichet, V. Legat, Adaptive high-order prediction of the drag correction factor for the upper-convected Maxwell fluid, J. Non-Newtonian Fluid Mech. 73 (1997) 95–114.