Lid-driven cavity flow of viscoelastic liquids

The lid-driven cavity flow is a well-known benchmark problem for the validation of new numerical methods and techniques. In experimental and numerical studies with viscoelastic fluids in such lid-driven flows, purely-elastic instabilities have been shown to appear even at very low Reynolds numbers. A finite-volume viscoelastic code, using the log-conformation formulation, is used in this work to probe the effect of viscoelasticity on the appearance of such instabilities in two-dimensional lid-driven cavities for a wide range of aspect ratios (0.125<height/length<4.0), at different Deborah numbers under creeping-flow conditions and to understand the effects of regularization of the lid velocity. The effect of the viscoelasticity on the steady-state results and on the critical conditions for the onset of the elastic instabilities are described and compared to experimental results.


Introduction
The fluid motion in a box induced by the translation of one wall -the so-called "lid-driven cavity flow" -is a classic problem in fluid mechanics [1] . The geometry is shown schematically in Fig. 1 and conventionally comprises a two-dimensional rectangular box of height H and width L of which the top (horizontal) wall -the "lid" -translates horizontally at a velocity U (for single phase fluids flowing isothermally the exact choice of moving wall is unimportant). For Newtonian fluids in rectangular boxes, the problem is governed by two dimensionless parameters: the Reynolds number Re ( ≡ρUH / η) and the aspect ratio ( ≡H / L ), where ρ is the fluid density and η is the dynamic viscosity. For very low Reynolds numbers -creeping-flow conditions, where Re → 0 -and low aspect ratios ( < 1.6 [2] ) a main recirculating region of fluid motion is induced which, due to the linearity of the Navier-Stokes equations under Stokes flow conditions, is symmetric about the vertical line x / L = 0.5, as shown in Fig. 1 . In addition to this main recirculation, smaller Moffatt [3] or corner eddies are also induced at the bottom corners (labelled "C" and "D" in Fig. 1 ): in fact at the bottom corners there is an infinite series of these vortices of diminishing size and strength as the corner is approached [3] . The * Corresponding author. Tel.: + 351 225081680. primary corner eddies grow in size with increasing aspect ratio and, at a critical aspect ratio of about 1.629 [2] , merge to form a secondary main cell albeit of much smaller intensity than the primary main cell (Ref. [1] gives a stream function decay ratio of 1/357 for these main cells). For higher aspect ratios this process repeats and additional main cells are created. Thus for "high" aspect ratios, essentially > 1.6, the main fluid motion near the translating wall (the main cell) is essentially unaffected by the aspect ratio. For a fixed aspect ratio -the majority of studies use = 1 -increasing the Reynolds number firstly breaks the foreaft symmetry about the vertical line x / L = 0.5 and then simulations reveal increasing complexity [1,4] . Since experiments show [5] that above a critical Reynolds number of around 500 the flow becomes three-dimensional and then time-dependent ( Re ≈ 825), we will not discuss these high Re steady two-dimensional simulations further here. Given the geometrical simplicity, combined with the rich observable fluid dynamics, the lid-driven cavity became a benchmark problem in the (Newtonian) fluid mechanics community for the development and validation of numerical schemes and discretization techniques [1,4,6,7] .
For viscoelastic fluids the literature is understandably less dense but, at least under creeping-flow conditions, much has been revealed by the limited number of studies to date.   [20] 1 .0 FV, 3D, Log conformation technique, CUBISTA Comminal et al. [21] 1 .0 Oldroyd-B, β = 0 .5 0.25 to 10 u ( x ) = 16 Ux 2 (1 − x ) 2 FD/FV, Logconformation, stream function Martins et al. [22] 1 .0 [23] 1 .0 Oldroyd-B, β = 0 .5  of the flow (which can be estimated as L / U ) and the Weissenberg number ( Wi ) which is defined as the ratio of elastic ( ∝ ληU 2 / H 2 ) to viscous stresses ( ∝ ηU / H ). Therefore De = λU / L and Wi = λU / H . Thus only in the unitary aspect ratio case are the two definitions identical: otherwise they are related through the aspect ratio ( De = Wi ). Experimentally, the papers of Pakdel and co-workers [8][9][10] detailed the creeping flow of two constant-viscosity elastic liquids (known as Boger fluids [11] ) -dilute solutions of a high molecular weight polyisobutylene polymer in a viscous polybutene oil -through a series of cavities of different aspect ratios (0.25 ≤ ≤ 4.0). Pakdel et al. [9] initially characterized the flow field at low wall velocities where the base flow remained steady and approximately two-dimensional: viscoelasticity was seen to break the fore-aft symmetry observed in Newtonian creeping flow and the eye of the recirculation region moved progressively further to the upper left quadrant (i.e. towards corner A of Fig. 1 ) of the cavities (incidentally, the fore-aft symmetry breaking due to inertia moves the eye towards corner B). At higher wall velocities, it was found [8,10,12] that the flow no longer remains steady but becomes time dependent. As inertia effects are vanishingly small in these highly-viscous Boger fluid flows, Pakdel and McKinley [10] associated this breakdown of the flow to a purely-elastic flow instability. The effect of cavity aspect ratio on this critical condition was found experimentally to occur at an approximately constant Deborah number (or, equivalently, 1/ Wi scaling linearly with aspect ratio ). Pakdel and McKinley [8] were able to explain this dependence on aspect ratio via a coupling between elasticity and streamline curvature and proposed a dimensionless criterion (now often referred to as the "Pakdel-McKinley" criterion [13][14][15] ) to capture this dependence on aspect ratio for both the lid driven cavity [8] and for a range of other flows [12] .
Grillet et al. [16] used a finite element technique to compute the effect of fluid elasticity on the flow kinematics and stress distribution in lid driven cavity flow, with a view to better understand the appearance of purely elastic instabilities in recirculating flows. In an effort to mimic the experiments and reduce, or circumvent, the numerical problems associated with the presence of a corner between a moving wall and a static one, the corner singularities were treated by incorporating a controlled amount of leakage. The results captured the experimentally observed upstream shift of the primary recirculation vortex and, concerning elastic instabilities, a dual instability mechanism was proposed, depending on aspect ratio (see also [17] ). For shallow aspect ratios, the downstream stress boundary layer is advected to the region of curvature at the bottom of the cavity resulting in a constant critical Weissenberg number; in deep cavities, the upstream stress boundary layer is advected to the region of curvature near the downstream corner (B in Fig. 1 ), resulting in a constant critical Deborah number.
Much as has been done for Newtonian fluids [6] , a number of studies have used the lid-driven cavity set-up to numerically test novel approaches or benchmark codes for viscoelastic fluids. These include Fattal and Kupferman [18] , who first proposed the log-conformation approach, and Pan et al. [19] , Habla et al. [20] , Comminal et al. [21] , and Martins et al. [22] , who later used that approach under the original or a modified form. Recently Dalal et al. [23] analysed the flow of shear-thinning viscoelastic fluids in rectangular lid-driven cavities, but also used the Oldroyd-B model for validation of the code. All these authors have applied a 4th order polynomial velocity regularization (what we will refer to later in Section 3 as "R1") to simulate the Oldroyd-B flow in 2D lid-driven cavities (Ref. [20] has tackled the 3D flow). In Table 1 we provide an overview of these previous numerical studies including the constant-viscosity viscoelastic model used (primarily Oldroyd-B with solvent-to-total viscosity ratio β = 0.5), the numerical and regularization methods used and the Weissenberg number reached. An exception was Yapici et al. [24] who solved the Oldroyd-B model for 0 ≤ Wi ≤ 1 with a finite-volume method, using the first-order upwind approximation for the viscoelastic stress fluxes in the rheological equation, and without recourse to the log-conformation approach. In marked contrast to all other studies with constant-viscosity viscoelastic models, including the present one, Yapici et al. [24] claim to be able to simulate viscoelastic lid-driven cavity flow without recourse to wall regularization.
Not surprisingly, there are a number of other studies with different types of non-Newtonian fluids, e.g. concerned with viscoplastic fluids where the interest is to identify the un-yielded central region (Mitsoulis and Zisis [25] , Zhang [26] , amongst others), and we shall use their results, in the Newtonian limit, as a basis for comparison.
In this work we re-visit the lid-driven cavity flow of viscoelastic fluids and investigate in detail how the choice of velocity regularization affects the viscoelastic simulations and the critical conditions under which the purely-elastic instability occurs.

Governing equations and numerical method
We are concerned with the isothermal, incompressible flow of a viscoelastic fluid flow, and hence the equations we need to solve are those of conservation of mass and of momentum together with an appropriate constitutive equation for the extrastress tensor, τ. In the current study we use both the upperconvected Maxwell (UCM) and Oldroyd-B models [27] τ = τ s + τ p , where λ is the fluid relaxation time, η s and η p are the solvent and polymer viscosities respectively, both of which are constant in these models (for the UCM model the solvent contribution is removed, η s = 0).
In order to increase the stability of the numerical method, we used the log-conformation procedure [28] , in which we solve for where A is the conformation tensor, which can be related to the extra-stress tensor as A = λ ηp τ p + I , E and R are the traceless extensional component and the pure rotational component of the veloc- , and I is the identity matrix [28,29] . An implicit finite-volume method was used to solve the governing equations. The method is based on a time-marching pressure-correction algorithm formulated with the collocated variable arrangement as originally described in Oliveira et al. [30] with subsequent improvements documented in Alves et al. [31] . The interested reader is referred to Afonso et al. [29] for more details and the corresponding numerical implementation and we only give a succinct overview here to avoid unnecessary repetition. The governing equations are transformed first to a generalized (usually non-orthogonal) coordinate system but the Cartesian velocity and stress components are retained. The equations are subsequently integrated in space over control volumes (cells with volume V P ) forming the computational mesh, and in time over a time step ( δt ), so that sets of linearized algebraic equations are obtained, having the general form: to be solved for the velocity components and for the logarithm of the conformation tensor. In these equations a F are coefficients accounting for advection and diffusion in the momentum equation and advection for the logarithm of the conformation tensor equations. The source term S φ is made up of all contributions that are not included in the terms with coefficients. The subscript P denotes the cell under consideration and subscript F its corresponding neighbouring cells. The central coefficient of the momentum equation, a P , is given by while for the log-conformation tensor equation is given by where a θ F contains only the convective fluxes multiplied by λ/ ρ.
After assembling all coefficients and source terms, the linear sets of Eq. ( 5 ) are solved initially for the logarithm of the conformation tensor and subsequently for the Cartesian velocity components. In general, these newly-computed velocity components do not satisfy continuity and therefore need to be corrected by an adjustment of the pressure differences which drive them. This adjustment is accomplished by means of a pressure-correction field obtained from a Poisson pressure equation according to the SIM-PLEC algorithm [32] to obtain a velocity field satisfying continuity. To discretize the convective fluxes, the method uses the CUBISTA high-resolution scheme, especially designed for differential constitutive equations [31] . In the current study our interest is restricted to creeping-flow conditions (i.e. Re → 0) in which case the advection term of the momentum equation (i.e. the second term on the left side of Eq. (2) ) is neglected. For the time-step discretization an implicit first-order Euler method was used, since we are primarily interested in the steady-state solution. We note that when steady-state conditions are achieved, the transient term used in the momentum equation (∂ u /∂ t ) for time marching vanishes, and we recover the Stokes flow equation for creeping flow.

Geometry, computational meshes and boundary conditions
The lid-driven cavity is shown schematically in Fig. 1 . To investigate the role of aspect ratio ( = H / L ) we have modelled eight different geometries ( = 0.125, 0.25, 0.5, 1.0, 1.5, 2.0, 3.0 and 4.0): thus "tall" enclosures correspond to aspect ratios greater than one and squat (or "shallow") enclosures to aspect ratios less than one. For each aspect ratio three consistently refined meshes have been used to enable the estimation of the numerical uncertainty of the results: for the square cavity additionally, a fourth finer mesh is used. For each geometry the central core of the cavity is covered with a uniform mesh which is progressively refined outside of this core region in both x and y directions so that the minimum cell size occurs in the four corners of the geometry. By construction, each mesh is symmetric about both the vertical and horizontal centrelines. Each mesh has an odd number of cells in both directions so that the variables are calculated exactly along the centrelines. The refinement procedure consists of halving the size of cells in both directions (and reducing the cell expansion/contraction factors accordingly) thus the total number of cells increases by essentially a factor of four between two meshes (to ensure an odd number of cells in each mesh the increase is not exactly a factor of four). The main characteristics of the meshes used for each aspect ratio are given in Table 2 , including the total number of cells (NC) and the minimum cell spacing which occurs at the corners ( x min / L and y min / L ).
The boundary conditions applied to the three stationary walls are no slip and impermeability (i.e. u = v = 0, where u and v are Table 2 Main characteristics of the computational meshes (NC = number of cells).
gives major numerical issues at the corners for viscoelastic fluids, due to the localized infinite acceleration applied to the fluid. The local extensional rate d u/ d x is infinite, theoretically, and the classical viscoelastic models here considered develop infinite stresses. Even though the numerical approximation introduces some degree of local smoothing, the method is unable to cope with the local stress peaks developed at the corners and fails to give a converged iterative solution. Hence, using this unregularized profile (which we shall refer to as "R0" henceforth) we could obtain converged solutions only in the case of low Weissenberg numbers (i.e. essentially Newtonian fluids only; for Wi > 0.02 there are already noticeable oscillations of the computed strength of the main recirculation). We note that this is in marked contrast with the results of Yapici et al. [24] who obtained results up to Wi = 1.0 for an unregularized profile. One common way of regularizing the lid-velocity is to use a polynomial function such that both the velocity and the velocity gradient vanish at the corners [19,20] . The use of such a regularization significantly reduces the strength of the main recirculation region within the cavity ( | ψ min | decreases in the Newtonian = 1 case by about 16% for example [33] ) and so to better mimic the unregularized idealized problem we also investigated the use of two weaker forms of regularization such that the velocity is uniform over the middle 60% of the moving wall R2 : and over 80% of its length R3 : Note that the velocity and velocity gradient also vanish at the corners for regularizations R2 and R3. However, although the velocity profile is continuous, the velocity gradient is not continuous at the points of change between the polynomial and the constant velocity profile, for R2 and R3. The different wall velocity profiles (i.e. R0, R1, R2 and R3), normalized using the peak velocity U , are shown in Fig. 2

Comparison with literature results and numerical accuracy
The numerical investigation of eight aspect ratios -using three different lid-velocity regularizations -for viscoelastic fluids over a range of Deborah (or Weissenberg) numbers, even in the limit of creeping flow, results in a large data set (approximately 600 simulations). Therefore, in this section only some representative data, which highlight typical levels of uncertainty, are presented.
Comparison of our data for Newtonian fluids with results in the literature (the square cavity case, R0), presented in Table 3 , shows excellent agreement. The minimum stream function value (i.e. the volumetric rate per unit depth of flow induced in the main recirculation region) agrees with values in the literature to within 0.05%. The minimum u velocity along x / L = 0.5 also agrees to literature results within 0.05%. There is a mild discrepancy ( ∼1%) with the    [4] in the minimum value of v computed along y / H = 0.5, but this may be a consequence of their "leaky" boundary conditions near the corners. The agreement of this quantity with the results of Yapici et al. [24] is, as for the other quantities, better than 0.1%. Comparison of our data for viscoelastic fluids, in this case for the Oldroyd-B model with a solvent-to-total viscosity ratio ( β = η s /( η s + η p )) of 0.5 and De = 0.5 and 1.0 (square cavity case R1), shown in Table 4 , indicates that, at De = 0.5, our results are in good agreement with Pan et al. [19] . At higher levels of elasticity, De = 1.0, the large normal stresses reveal a greater degree of sensitivity and non-negligible differences between results in both our finest two meshes and, also, in comparison with the data of Pan et al. [19] . The effects of mesh density on the accuracy of the numerical results are shown in Tables 5-7 for = 1, 0.125 and 4, respectively. Overall, the difference between the results on the finest mesh and the extrapolated values -obtained using Richardson's technique [34] -are less than 0.2% for these quantities.

Creeping Newtonian flow
The streamline patterns for Newtonian flow -wall regularization R3 and mesh M3 -are shown in Fig. 3 for low aspect ratio ( < 1, Fig. 3 (a)) and high aspect ratio ( ≥ 1 Fig. 3 (b)). As discussed in the Introduction (and in [2] ), for high aspect ratios the streamlines essentially collapse in the top region of the cavities and the maximum absolute value of the stream function -the variation of which with aspect ratio is shown in Fig. 4 -becomes independent of aspect ratio for ≥ 1. For small aspect ratios the scenario is more complex although the asymptotic limit when → 0 does allow analytical expressions for the maximum stream function, velocity and stress components to be derived (presented Table 5 Effect  in Appendix A ). Under these assumptions the maximum absolute value of stream function is expected to vary linearly with the aspect ratio as | ψ min | UL = 4 27 , and this linear relationship is also included in Fig. 4 . Excellent agreement between this analytical solution and computations can be seen with aspect ratios up to = 0.5, especially for the unregularized and weakly regularized   Reg wall velocities (R0 and R3). In addition, the intersection between that linear variation and the value | ψ min | /UL = 0 . 1 at high aspect ratio shows that the "tall" cavities start at ≈ 0.7. The effects of wall regularization are subtle, yet important. In the square cavity case, = 1, regularizing using the standard polynomial function [18,19] (R1 in our nomenclature) reduces quantitatively the strength of the main recirculation, by about 16%, (in agreement with previous studies [33] ) but the qualitative effect on the streamlines ( Fig. 5 ) appears to be minor except close to corners A and B, as might be expected. At the lower aspect ratios, however, there are significant qualitative differences between the streamline patterns and the unregularized streamlines are essentially straight over the middle 75% of the cavity. Changing the regularization such that it better approximates the unregularized case, e.g. R3, essentially increases the vortex strength back to its unregularized value (see Table 5 for example) and better captures the streamline patterns (although close to the corners A and B differences are still apparent - Fig. 5 ). To better illustrate these effects, in Fig. 6 we plot contours of the flow type classifier. The flow-type parameter ξ is used to classify the flow locally, and here we use the criterion proposed by Lee et al. [35] , ξ ≡ | D |−| | | D | + | | , where | D | and | | represent the magnitudes of the rate of deformation tensor and  vorticity tensor, | D | = 1 2 ( D : D ) and | | = 1 2 ( : T ) . As such, ξ = 1 corresponds to pure extensional flow, ξ = 0 corresponds to pure shear flow and ξ = −1 corresponds to solid-body rotation flow.
As seen in Fig. 6 , the regularization of the lid velocity significantly affects the flow close to the lid ends. Without regularization (R0), the flow close to the wall is mainly shear dominated. However, the wall regularization induces a strong component of extensional flow close to corners A and B due to the acceleration and deceleration of the fluid at these corners. As this acceleration region decreases (R1 → R2 → R3) the region of extensional-dominated flow also decreases and approaches the "true" lid-driven cavity flow field (R0).
Given the basic modification to the flow field induced by the regularization classically used for viscoelastic fluids [18,19] , care must be taken when comparing regularized simulation results with experimental results [9,10,17] . This issue is probably why Grillet et al. [16] implemented a "leaky" boundary condition at corners A and B. In contrast here we attempt to tackle this problem via modification of the classical regularization (R1).

Steady-state flow field
The addition of fluid elasticity induces changes to the flow pattern in the lid-driven cavity, and, in particular, a breaking of fore-aft symmetry relative to the x / L = 0.5 line. The large normal stresses that are generated for the viscoelastic fluid as De increases are advected in the downstream direction leading to an increase of the flow resistance, and to compensate for this effect the eye of the recirculation region progressively shifts in the upwind direction -towards corner A as illustrated in Fig. 1 breaking the symmetry observed for Newtonian fluids. This effect is in excellent agreement with experimental observations for Boger fluids [9] . The increase of the normal stresses with increasing De , and concomitant higher flow resistance, induces a decrease in the strength of the main recirculating flow (the flow closest to the lid), i.e. a reduction of | ψ min | as illustrated in Fig. 8 for different aspect ratios as a function of De (for three lid velocity regularizations). This effect is akin to the vortex suppression by elasticity seen in many other flow situations, for example in sudden expansion geometries, and entails the coupling of elastic hoop stresses and curved streamlines (as discussed below in Section 6.2 ). For the lower aspect ratio cases and the regularizations (R2 and R3) closer to the unregularized situation (R0), the streamlines in a large region about the central section of the cavity are straight, with curvature confined to a small region near the lateral walls (see e.g. Fig. 5 top), so the mechanism of the vortex suppression should be less effective. Indeed, we find in Fig. 8 a slight initial increase, although small (about 1-2%) of the vortex intensity with elasticity for the cases of = 0.125 and 0.25 with regularizations R2 and R3, which we interpret as resulting from straighter streamlines.
Similarly to the Newtonian fluid flow case, the intensity of recirculating fluid increases with aspect ratio up to = 1, before saturating, and further increases of have a negligible effect on | ψ min | as shown by the data collapse for ≥ 1. The additional recirculation zones that are formed as increases, e.g. shown for a Newtonian fluid in Fig. 5 , are significantly weaker, having a negligible effect on | ψ min | . With the regularization R1 (black symbols), since the lid-velocity profile is smoother and has a lower average velocity relative to the R3 regularization, the maximum value of the stream function magnitude is also lower and higher Deborah numbers can be achieved prior to the onset of a purely-elastic instability, which we discuss next.

Onset and scaling of a purely-elastic instability
The critical conditions for the onset of a purely-elastic instability are presented in Fig. 9 , both in terms of a critical reciprocal Weissenberg number (1/ Wi cr ) and a Deborah number ( De cr ) as a function of aspect ratio, using three different lid regularizations (R1-R3), computed using mesh M3. The critical condition is identified when a steady-state solution can no longer be obtained: thus the purely-elastic instability, in all cases examined here, gives rise to a time varying flow in agreement with experimental observations [10] . We confirmed that this critical condition is independent of the time-step used in the time-marching algorithm.
For a given level of wall regularization, the flow instabilities for different aspect ratios occur approximately at a constant Deborah number, e.g. for R1 the values are within De = 0.625 ± 0.050. Consequently, 1/ Wi cr varies linearly with aspect ratio, with 1/ Wi cr = 1.60 for regularization R1. As the regularization is weakened, and the forcing better approximates the "true" unregularized case i.e. R1 → R2 → R3, the critical De (or Wi ) decreases sig- nificantly from 0.63 (R1) to 0.33 (R2) to 0.18 (R3). Thus the precise regularization used can decrease the critical De to about one third. The use of an average wall velocity ( U = L 0 u dx /L ), instead of the peak velocity U , reduces these differences in critical value slightly (to about half). The use of a more appropriate characteristic time, based on the distance over which the lid velocity grows, to define De * and Wi * as described in Section 3 , is much better able to collapse the critical conditions as shown in Fig. 9 (b). For the two weakest forms of regularization (R2 and R3) this practically collapses the critical values. If this scaling is representative of the controlling dynamics it would imply that the unregularized lid is unstable for vanishingly small elasticities (given the infinite acceleration, or zero acceleration time, for a fluid element to go from rest to velocity U, or equivalently De * → ∞ ), but as shown in Section 3 , numerically De * is finite since the mesh elements do not have zero length. Fig. 10 (a) presents the contours of the Pakdel-McKinley criterion M crit = ( λu/ ) ( τ 11 / η 0 ˙ γ ) [8,12] , (where τ 11 is the tensile stress in the local streamwise direction, ˙ γ is the local shear rate and is the local streamline radius of curvature), for = 0.25, 0.5, 1, 2 and 4 for the highest steady De . The maximum values of M crit are of the same order as observed in previous studies, between 3 and 4 [12] , namely 3.6 ( = 0.25), 3.9 ( = 0.5), 3.5 ( = 1, 2 and 4). These critical M values are located near the downstream corner, where large normal stresses generated on approaching corner B are advected into a region of high streamline curvature.
Finally, contours of the flow type parameter are presented in Fig. 10 (b), close to the highest stable De for which the flow remains steady for = 0.125 and 1, for both regularizations R1 and R3. For = 1, the flow is mainly rotational close to the main vortex centre, shear dominated close to the lid and highly extensional near the top corners in a direction at 45 º and in the lower part of the cavity, albeit here the deformation rates and hence stresses are always modest. For = 0.125, the flow is mainly shear dominated close to the lid and bottom wall and highly extensional close to the corners and along a thin strand at y / H ≈ 1/3.

Conclusions
Viscoelastic creeping flow in a lid-driven cavity was analysed numerically using a finite-volume numerical methodology in combination with the log-conformation technique. The effects of aspect ratio, strength of elasticity via the Deborah/Weissenberg numbers and, in contrast to previous studies, the effect of regularization type for the lid velocity were investigated. As discussed in the introduction and as observed in previous studies for Newtonian fluids, the streamlines essentially collapse and the stream function becomes independent of aspect ratio, for ≥ 1. For low aspect ratios, the maximum value of the stream function magnitude agrees with the analytical solution for a simple parallel-flow approximation.
The effect of elasticity on the steady flow characteristics was elucidated, including the experimentally-observed shift of the main vortex centre in the direction of the upstream corner (A) and the breaking of fore-aft symmetry with increasing elasticity. Increasing the elasticity was also found to reduce the flow strength induced by the lid. The critical conditions for the onset of purely-elastic flow instabilities were characterized by plotting 1/ Wi cr (and De cr ) as a function of aspect ratio. The value of the Pakdel-McKinley criterion [12] at the critical conditions was similar to previous studies, between 3 and 4. In accordance with the experimental results of Pakdel and McKinley [8,10,12] , we find a linear scaling of reciprocal Weissenberg number with aspect ratio, corresponding to a constant Deborah number, the numerical value of which is dependent on the wall regularization used. By introducing a modified Deborah number, based on an average acceleration time for the flow adjacent to the lid to accelerate from u = 0 to u = U , a reasonable collapse for the various is achieved. The flow-type parameter was also computed to characterize the different regions of the lid-driven cavity flow. Further studies are required to investigate more closely the generated time-dependent flows.
We can also use the velocity distribution to estimate the maximum dimensionless shear stress on the lid ( = η 0 ( d u /d y) ( y = H ) ) and also the maximum dimensionless axial normal stress