Finite volume simulation of viscoelastic laminar flow in a lid-driven cavity

J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

Contents lists available at ScienceDirect

Journal of Non-Newtonian Fluid Mechanics journal homepage: www.elsevier.com/locate/jnnfm

Finite volume simulation of viscoelastic laminar flow in a lid-driven cavity Kerim Yapici a , Bulent Karasozen b , Yusuf Uludag a,∗ a b

Department of Chemical Engineering, Middle East Technical University, Inonu Bulvari, 06531 Ankara, Turkey Department of Mathematics and Institute of Applied Mathematics, Middle East Technical University, 06531 Ankara, Turkey

a r t i c l e

i n f o

Article history: Received 3 November 2008 Received in revised form 1 June 2009 Accepted 5 August 2009

Keywords: Lid-driven cavity Oldroyd-B fluid Finite volume Collocated grid

a b s t r a c t A finite volume technique is presented for the numerical solution of steady laminar flow of Oldroyd-B fluid in a lid-driven square cavity over a wide range of Reynolds and Weissenberg numbers. Second order central difference (CD) scheme is used for the convection part of the momentum equation while first order upwind approximation is employed to handle viscoelastic stresses. A non-uniform collocated grid system is used. Momentum interpolation method (MIM) is used to evaluate face velocity. Coupled mass and momentum conservation equations are solved through iterative SIMPLE (semi-implicit method for pressure-linked equation) algorithm. Detailed investigation of the flow field is carried out in terms of velocity and stress fields. Differences between the behavior of Newtonian and viscoelastic fluids, such as the normal stress effects and secondary eddy formations, are highlighted. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Most materials used in many industries such as plastic, food, pharmaceuticals, electronics, dye, etc. exhibit viscoelastic properties under their processing or flow conditions. Viscoelasticity stems from partial memory of the material on the deformation history. The memory fades away over material dependent relaxation time, which is at least comparable to the observed time scale in the flow. Due to the elasticity of such materials, their hydrodynamic behavior differs from simple Newtonian fluids in many important respects. Rod climbing, siphoning, secondary flows are all common examples to how a viscoelastic fluid can exhibit quite different flow behavior than a Newtonian fluid would do under similar flow conditions. In industrial processes involving flow of viscoelastic materials, understanding complexities associated with viscoelasticity can lead to both design and development of hydrodynamically efficient processes and improved quality control of the final products. Number of experimental techniques has been employed to study influence of the viscoelasticity on the flow kinematics in complicated geometries [1–4]. With the advances in digital computer technology and with the numerical methods of enhanced accuracy, stability, efficiency and robustness, number of studies based on the computer simulation of the viscoelastic flows has steadily increased. Therefore, it has been possible to attack complex viscoelastic flows using computer simulations instead of experiments which can be labor intensive,

∗ Corresponding author. Tel.: +90 312 210 4374; fax: +90 312 210 2600. E-mail address: [email protected] (Y. Uludag). 0377-0257/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2009.08.004

expensive, hazardous and even impossible to conduct in a laboratory [5]. Computational analysis of a viscoelastic flow requires rheological models or constitutive equations that provide relation between stress and material deformation rates [6]. Many constitutive equations have been developed to describe the behavior of viscoelastic fluids. Some of the non-linear differential type equations are Upper Convected Maxwell (UCM), Oldroyd-B, Phan-Thien Tanner (PTT) and White-Metzner models. These are commonly used in the numerical simulations. So far many researches have been conducted to develop robust and stable numerical algorithm to solve viscoelastic fluid flows at moderately high values of Weissenberg number, We, which has been one of the main challenges in the simulations. The loss of convergence of numerical simulation at some critical value of Weissenberg number is called in literature as high Weissenberg number problem [7–9]. It is known in the literature that at high We most of the numerical methods fail to converge [8,10,11]. Common causes are: extremely large stress components associated with abrupt change of boundary [12,13], use of unsuitable boundary conditions [11], and changing equation characteristic from hyperbolic to elliptic or hyperbolic to parabolic when inertial terms are considerable [14,15], and constitutive models [8]. So far various numerical techniques have been developed to predict flow kinematics of viscoelastic fluids accurately by using number of rheological models. Finite difference [16–18], finite element [19,20], spectral finite element [21] and finite volume [6,22–26] techniques have been used in the discretization of the non-linear flow systems. Finite difference method applications are generally restricted to simple geometries. Finite volume approximation of fluid flow systems, on the other hand, can be

52

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

advantageous in terms of computer space and time requirements as well as in terms of numerical stability compared to the finite element method [26]. Source term in the momentum equation includes non-linearly coupled extra stress terms. In the discretization step, these terms can be handled by coupled or decoupled strategy. Description of both approaches is reported by Tanner and Xue [27]. Moreover implementation of decoupled strategy for finite volume method in steady and transient discretization form can be found elsewhere [23,27]. Most numerical simulations of viscoelastic fluids using finite volume technique are based on staggered grids arrangements [6,23,24,28]. Here pressure and normal stress components  xx and  yy are located at the center of the control volume, while generally  xy is located at corners and velocities are placed on the faces of the control volume. To avoid checkerboard or zig-zag pressure distribution, staggered grid arrangements are also preferred in Newtonian flow simulations. The other kind of arrangement is non-staggered grid in which all flow variables are located at the center of the control volume. These two kinds of grid arrangements are compared in detail by Peric et al. [29] by performing three test cases involving a Newtonian fluid; lid-driven cavity, backward facing step and flow through a pipe with sudden contraction. Some advantages of collocated grids are listed as follows [29]: (a) because all variables stored at the same location, the coefficients in the discretization equations are identical for all velocities. The stress coefficients in the equations are also the same. This feature reduces computational time and required memory when higher order schemes like quadratic upstream interpolation for convective kinematics (QUICK) discretization method are employed in both Newtonian and viscoelastic flow simulations. (b) Treatment of the boundary conditions is also easy. (c) Implementation of the multigrid method is simpler than that of staggered since the same interpolation factor is used for all variables. However, non-staggered grids produce checkerboard pressure with coarse grids [30]. To overcome this problem, Rhie and Chow [31] first proposed the well-known and popular momentum interpolation method (MIM). Although this approach removes the undesired pressure distribution, under-relaxation parameter dependency occurs [32]. Majumdar [32] modified the MIM of Rhie and Chow [31] to get around this problem. Detailed derivations of the proposed algorithm of Ref. [32] along with two new momentum interpolation methods are well documented by Yu et al. [33]. Simulations based on finite volume method with non-staggered grids are also used in the simulations of flows of viscoelastic fluids [25,34–36]. Oliveira et al. [25] performed MIM interpolation method to remove pressure–velocity–stress decoupling for orthogonal and non-orthogonal graded meshes. They employed finite volume method with SIMPLEC algorithm for slip-stick and flow around a circular cylinder cases with a UCM fluid. Several schemes have been developed and used to approximate the convection term, such as hybrid differencing scheme [37], power-law differencing scheme [38], skew upwind differencing scheme (SUDS) [39], central difference (CD) scheme [40], quadratic upstream interpolation for convective kinematics (QUICK) [41] and sharp and monotonic algorithm for realistic transport (SMART) [42]. Most of the computational studies involving viscoelastic materials consider sudden expansion or contraction flows. Few studies on the lid-driven cavity geometry exist in the literature [43–50]. Despite its simple geometry, lid-driven cavity offers many challenges that are encountered in industrial applications in the form of singularity points, circulations and bifurcations in the flow. Therefore it is highly likely that a numerical methodology efficiently simulating flow of viscoelastic material in this benchmark geometry, may handle complex flows associated with the industrial applications such as agitation.

Grillet et al. [47] used finite element method to analyze stability of the recirculating flows in two-dimensional lid-driven cavity geometry. In order to eliminate the effect of the corner singularities they introduced small leakages at the upstream and downstream corners of the cavity. Setting the amount of leakage was critical due to its strong impact on the stresses and kinematics at the lid corners hence on the entire flow field. Setting the amount of leakage, however, suffers from experimental verification making the value used in the study ambiguous. They also investigated the effect of the cavity aspect ratio on the flow field and proposed dual elastic instability criterion. They suggested a mechanism for elastic instability based on the convective elastic stresses either from upstream corner to the down stream corner or from downstream corner to the lid bottom depending on the height/length ratio of the cavity. Fattal and Kupferman [49] employed second order finite difference scheme and converted the algebraic equations into logarithmic form which enabled them to decrease the impact of the Weissenberg instability. They utilized this method for lid-driven cavity at creeping flow condition by using Oldroyd-B constitutive model. Oliveira [50] used finite volume based collocated mesh arrangements for the simulation of creeping flow of a FENE-CR fluid in a lid-driven cavity. To remove corner singularities effect, the lid velocity was expressed by describing polynomial velocity distribution which goes to zero at the corners. The FENE-CR fluid had solvent viscosity ratio of 0.79. He was able to get converged solutions up to Deborah numbers of 10. Pak et al. [4] investigated both laminar and turbulent flow of viscoelastic fluid through a circular pipe with a sudden expansion. They observed that reattachment length for viscoelastic fluids in the laminar regime was shorter than that of Newtonian fluid. On the other hand, in turbulent flow regime, the reattachment length for viscoelastic fluids turned out to be two or three times longer. The hydrodynamic behavior of viscoelastic fluids in lid-driven cavity was also reported in the literature for creeping flow cases. However, the effect of the Reynolds number on the flow field in the literature is not available to the best knowledge of the authors. In the present study, the influence of the fluid elasticity is examined computationally in a square lid-driven cavity with and without inertial effects. In our numerical analysis, viscoelastic behavior is modeled by the Oldroyd-B constitutive equation which predicts a quadratic first normal stress difference and constant viscosity. No artificial leakage, hence uncertainties, are introduced. The simulated flow system with boundary condition is depicted in Fig. 1. The convective terms in the momentum equations are discretized by second order central differences. While upwind is used in the constitutive equations for the stresses the semi-implicit method for the pressure-linked equation (SIMPLE) is used to solve the coupled continuity, momentum and constitutive equations. This paper is organized as follows. Description of the governing equations is presented in Section 2. In the subsequent section, used numerical tools are outlined. In Section 4, code validation and analysis of the numerical results that are obtained for lid-driven cavity are included along with comparison with Newtonian fluid results. In the last section, conclusions and recommendations are presented.

2. Governing equations We consider steady, incompressible and isothermal flow of Oldroyd-B [51,52] fluid in two-dimensional Cartesian coordinate system (x,y). The followings are dimensionless form of the continuity, momentum and constitutive equations used in this study. Detail of the derivations associated with the constitutive equation

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

53

Fig. 2. Location of the flow variables on non-uniform collocated grids. Fig. 1. Square lid-driven cavity geometry with boundary conditions.

Set of algebraic equations which are non-linear due to the source term in the constitutive equations. To make the equations linear, first, source term S is assumed to be a linear function of variable

can be found elsewhere [51,52]. ∂u ∂v + =0 ∂x  ∂y  ∂ ∂u ∂ Re uu − (1 − wr) + ∂x ∂x ∂y ∂ ∂x



Re uv − (1 − wr)

∂v ∂x



+

∂ ∂y

 Re vu − (1 − wr)

 Re vv − (1 − wr)

∂u ∂y

∂v ∂y

 =−

 =−

∂xy ∂xx ∂p + + ∂y ∂x ∂x

∂yy ∂xy ∂p + + ∂y ∂y ∂x

Dimensionless form of the stress components  xx ,  yy and  xy are given as follows:



∂v ∂u − ∂y ∂x



xx +

∂ ∂ (We uxx ) + (We vxx ) = We ∂x ∂y

yy +

∂ ∂ (We uyy ) + (We vyy ) = We ∂x ∂y

xy +

∂ ∂ 1 (We uxy ) + (We vxy ) = − We(xx − yy ) 2 ∂x ∂y



∂v ∂u − ∂x ∂y

xy + 2We

∂u xx + We ∂x

xy + 2We

∂v yy + We ∂y

 

∂v ∂u − ∂y ∂x





+ wr

where wr = 1 − ˇ and the parameter ˇ is the ratio of the retardation and relaxation time and is defined by ˇ = 2 /1 . The Reynolds number and the Weissenberg number which is defined as the ratio of characteristic fluid relaxation time to characteristic time scale in the flow are given through Re = UH/ and We = 1 U/H respectively. 3. Numerical method In this study finite volume formulation [38,53] is used to solve conservation equations of mass, momentum, and constitutive equations that provide relation between stress and material deformation rates [6]. Continuity, momentum and constitutive equations can be written in general form as follows: ∂ ∂x



u −

∂ ∂x



+

∂ ∂y



u −

∂ ∂y



= S

(1)

(3)

where  is either density  or relaxation time , depending on if conservation or constitutive equation is considered;  is one of the dependent variables; is the diffusion coefficient and S is the source term.

 

∂v ∂u + ∂y ∂x ∂u ∂v + ∂y ∂x

∂v ∂u + ∂y ∂x

 xy + 2wr

∂u ∂x

xy + 2wr

∂v ∂x





+

1 We(xx + yy ) 2

(2)



∂v ∂u + ∂y ∂x



 such that, S = SC + SP P

(4)

where SC is constant part of the S that is independent of  while SP is the coefficient of P which is set as negative to enhance the numerical stability [38]. In our calculations, central difference is employed for both the approximation of the gradients which are diffusion terms in Eq. (3) and for the discretization of convective terms in the momentum equations. Upwind differencing scheme is used in the constitutive equations for the stresses. We used non-uniform collocated grid arrangements to discretize the governing sets of equations in the computing domain. In this method, all flow variables are located at the center of the control volumes as shown in Fig. 2. The final form of the two-dimensional discretized governing equation over the control volume (see Fig. 3) can be expressed symbolically as follows: AP P = AE E + AW W + AN N + AS S + b

(5)

54

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

Fig. 4. Planar channel flow geometry.

Fig. 3. Schematic diagram of a control volume.

Discretized form of the x-momentum and  xx equations are given as AuP ui,j = AE ui+1,j + AW ui−1,j + AN ui,j+1 + AS ui,j−1 + bu APxx xx i,j = AE xx i+1,j + AW xx i−1,j + AN xx i,j+1 +AS xx i,j−1 + bxx

(6)

Furthermore, gradients of the pressure and stresses in the momentum equations are obtained by linear interpolation. While gradients of velocities in constitutive equation are computed by central differences at interior domain and near the boundary they are calculated by introducing a polynomial function such as:

v=

ax2

+ bx + c

u = ay2 + by + c

∂v = 2ax + b ∂x ∂u = 2ay + b ∂y

technique is adopted in such a way that source term of the momentum equations is treated as pseudo-body forces that are already determined in the previous iteration [23]. The set of linearized algebraic equations are then solved using the Line-by-Line technique based on the TDMA (Thomas algorithm or the tridiagonal matrix algorithm) and the alternative direction implicit (ADI) scheme. For pressure-correction equation, four ADI sweeps are made [29] before correcting pressure and velocities. To stabilize the calculations, global under-relaxation factors are used depending on values of the Reynolds and Weissenberg numbers. For instance, at Re = 0, the factors are set as 0.4 and 0.7 for velocities and pressure, respectively. For the stress components the under-relaxation factor of 0.1 is employed for We ≤ 0.5. At larger We, it is reduced to 0.05 for all meshes. The solution process is reiterated until the maximum relative change of flow variables (u, v, p,  xx ,  yy ,  xy ) are less then a prescribed tolerance or residual as:

 res = MAX

|n+1 − n | |n+1 |

 ≤ 1 × 10−5

(8)

(7)

Discretized two-dimensional momentum and constitutive equations are written for interior nodal points. Boundary nodes, however, require special treatment [53]. In this study to eliminate checkerboard pressure distribution, cell face velocities are evaluated by momentum interpolation method (MIM) suggested by Majumdar [32]. MIM was first proposed by Rhie and Chow [31]. Detailed procedure of the original MIM and two new interpolation methods, called MMIM1 and MMIM2 were well documented by Yu et al. [33]. Yu et al. [33] recommended that MIM should be used not only in the mass-source evaluation but also for the evaluation of coefficient of the discretization equation. Hence, the modified MIM by Majumdar [32] is used to evaluate coefficients in the constitutive equations. Due to the collocated arrangement, velocity gradients in the source term of constitutive equation should be evaluated at the center of the control volumes. This is done by using cell face velocities that are evaluated using MIM, instead of velocities at nodal points. 3.1. Solution algorithm To date various algorithms have been developed to calculate pressure field. Semi-implicit method for the pressure-linked equation (SIMPLE) was the first proposed by Patankar and Spalding [54] and it has been widely used in the literature. After that several modifications have been made to enhance its efficiency, such as SIMPLE revised (SIMPLER) [38], SIMPLE consistent (SIMPLEC) [55] and pressure-implicit with splitting of operators (PISO) [56]. In this study, the SIMPLE is employed to solve the coupled system of the continuity, momentum and constitutive equations. Source terms of the momentum equations contain extra stresses which are nonlinearly coupled with the stress equations. Here, the decoupling

where  = (u, v, p,  xx ,  yy ,  xy )T .

4. Results and discussion All the simulations presented were carried out on single computer platform, including an AMD Athlon 64 FX-57 processor. The computer code, which we name as NONSOL (NOn-Newtonian flow SOLver), was written based on the aforementioned algorithms. The code is first used for developing flow of an Oldroyd-B fluid in a simple geometry, a planar channel as shown in Fig. 4, for which analytical solution can easily be obtained. Following satisfactory agreement between the simulation results and analytical solution for this geometry, steady lid-driven flow in a square cavity is considered. Details and results of these two different flow geometries are given in the following sections.

4.1. Flow in entrance of planar channel As shown in Fig. 4, we first consider steady creeping flow of an Oldroyd-B fluid model through a planar channel, which has length of 8H and width 2H (H = 1). At the inlet the flow is assumed to be plug flow except in the immediate vicinity of the wall, as described by Gaidos and Darby [20] and reformulated by Na and Yoo [52].

 u(y) =

1.0448 (−46.92 + 110.09y − 63.17y2 )

for y ≤ 0.8714 for y > 0.8714

(9)

Because of the hyperbolicity of the constitutive equations the stress components must be specified at the entrance [20]. They are calculated by using the assumed entrance flow distribution and by the

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

55

Fig. 5. Velocity profiles (a) at near the exit of the channel and (b) at centerline for Re = 0.

following relations.

 xx = 2We(1 − ˇ)

∂u ∂y

2 (10)

yy = 0 xy

∂u = (1 − ˇ) ∂y

No slip boundary conditions are assumed at the solid wall for velocity, u = 0,

∂v = 0, ∂x

v = 0,

∂v = 0, ∂y

∂u =0 ∂x

(11)

The stress components at the wall are obtained by solving constitutive equations and by applying known velocity boundary conditions [57] in Eq. (11). At the centerline ∂u = 0, ∂y

v = 0,

4.2. Lid-driven cavity flow

∂v =0 ∂x

(12)

Similarly the stress components at the symmetry line are obtained by solving constitutive equations and by using Eq. (12). At the exit the Neumann boundary conditions are imposed for the flow variables ∂u = 0, ∂x

∂v = 0, ∂x

∂xx = 0, ∂x

of Neumann type, Dirichlet boundary conditions are used for pressure-correction equations. We used uniform collocated grids of 65 points in y-direction. We have chosen polymer contributed viscosity as ˇ = 0.2 in the Oldroyd-B model. Flow conditions are set as Re = 1 × 10−4 after that this value is stated as 0 and We = 0.1, 0.5, 0.8. Numerical results of the developing velocity and stresses are shown in Figs. 5 and 6 along with their corresponding analytical solutions. Fig. 5 illustrates influence of the We on the fully developed (x > 7) velocity profiles and centerline velocity profiles at Re of 0. The maximum value of the velocity is 1.5 for all We numbers. These results also confirm that Oldroyd-B model does not exhibit shear thinning behavior. Profiles of the corresponding normal and shear stress are presented in Fig. 6. Numerical results show that excellent agreement between the simulation results and those of the analytical solutions is obtained.

∂yy = 0, ∂x

∂xy =0 ∂x

(13)

Boundary values for pressure-correction equations are determined by considering velocity boundary conditions [58]. For instance, if the velocity boundary conditions are of Dirichlet type, Neumann type boundary conditions are applied to the pressure-correction equations. Conversely, if the velocity boundary conditions are

In this section first, the lowest possible mesh resolution required to get mesh independent solutions is investigated. For that purpose stream function value, first normal stress difference and horizontal velocity at the center of cavity are considered at We = 0.1 and 1 at Re = 0. Then determined mesh size is employed in the rest of the calculations. Numerical results of the Oldroyd-B fluid having polymer viscosity, ˇ, of 0.3 are presented within We range from 0 to 1. 4.2.1. Boundary conditions for stresses Boundary values of stress components can be determined by solving the constitutive equations with the available velocity field.

Fig. 6. Comparison of computed and analytical profiles of (a) normal and (b) shear stresses at near the exit of the channel for Re = 0.

56

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

Table 1 Main characteristics of meshes used. Meshes

(2) Bottom wall boundary conditions parallel to the x-axis: In this case the conditions are as follows: Minimum x

Number of nodes 129 × 129 193 × 193 257 × 257 305 × 305

M1 M2 M3 M4

0.00524 0.00349 0.00262 0.00220

u = 0,

u = 1,

v = 0,

∂v = 0, ∂x

∂v =0 ∂y

 (14)

If we apply these conditions in the constitutive Eq. (2) we obtain the following non-linear equations:





∂ ∂u xy (We uxx ) = 2We ∂x ∂y ∂ (We uyy ) = 0 yy + ∂x     ∂ ∂u ∂u xy + (We uxy ) = Weyy + wr ∂x ∂y ∂y

∂v = 0, ∂y

∂u =0 ∂x

(16)

Stress boundary conditions for the bottom wall can then be expressed by the following equations.

(1) Top wall boundary conditions parallel to the x-axis: ∂u = 0, ∂x

∂v = 0, ∂x

v = 0,

xx = 2We wr yy = 0 xy = wr



∂u ∂y

∂u ∂y

2 (17)



xx +

(15)

Therefore, the values of  xx ,  yy and  xy on the lid are obtained by solving the set of equations above.

(3) Solid walls boundary conditions parallel to the y-axis (left and right wall): The values of  xx ,  yy and  xy at the solid walls are evaluated through the following conditions: u = 0,

v = 0,

∂u = 0, ∂y

∂v = 0, ∂y

∂u =0 ∂x

(18)

Fig. 7. Effect of the total number of nodes on the stream function value ( ), horizontal velocity (u) and first normal stress differences (N1 ) at location of x = 0.5 and y = 0.5: (a) We = 0.1 and (b) We = 1 at Re = 0.

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

57

Fig. 8. The stream functions for different We at Re = 0. Contour levels shown for each plot are the same.

xx = 0



yy = 2We wr



xy = wr

∂v ∂x



∂v ∂x

All the results in the following sections are obtained by using dense mesh in M4.

2 (19)

The boundary conditions on pressure-correction are of Neumann type for all boundaries.

4.2.2. Effect of mesh size on the flow variables Mesh refinement is one of the most important issues in viscoelastic flow simulations in order to acquire accurate results and capture realistic flow behavior near the singularities of lid-driven cavity or contraction flow problems. For that reason we examined effect of mesh refinement on convergence of the solutions by using various meshes sizes starting from coarse mesh 65 × 65 to dense mesh 305 × 305 with the increments of 16 in i and j directions. The main characteristics of selected meshes used in our simulations are tabulated in Table 1. In order to demonstrate mesh independence of the solutions, stream function value, first normal stress difference and horizontal velocity taken at a fixed location in cavity is shown in Fig. 7 as a function of total number of nodes at x = 0.5, y = 0.5 or at the center of the cavity. It shows that as number of nodes increases the flow variables tend to be mesh independent. Also depicted in the figure, mesh refinement study for viscoelastic fluid simulations should be done by considering at least three or more flow variables in order to ensure mesh independent results. This is especially highlighted in Fig. 7(a).

4.2.3. Results at creeping flow conditions At creeping flow conditions, the stream functions and vorticity at a fixed polymer viscosity, ˇ = 0.3 and different We are depicted in Figs. 8 and 9, respectively. We = 0 results are used as Newtonian solution in the text. Stream function, , and vorticity, ω, are obtained by solving the following equations after convergence is attained.

∇2

= −ω

ω=−

∂v ∂u + ∂y ∂x

(20) (21)

In Newtonian cavity simulations, stream function and vorticity values and their locations are usually taken as a code validation. At We = 0 and Re = 0, i.e. Stokes flow, fore-aft symmetry is observed, a typical Newtonian fluid behavior. However, for a viscoelastic fluid, as the value of We increases the primary vortex shifts in the upstream direction (see Fig. 10) and streamlines near the downstream corners become more curved due to the increased normal stress effects. This phenomenon was also observed experimentally by Pakdel et al. [45] using laser Doppler velocimetry (LDV) and digital particle image velocimetry (DPIV). The effect of the elasticity is more pronounced at the downstream than the upstream due to the finite relaxation of the normal stresses as shown by the vorticity contours in Fig. 9. Fig. 11 shows the effect of We on the horizontal (u) and vertical (v) velocity components distribution at x = 0.5 and y = 0.5, respectively. With increasing We, maximum magnitudes of the velocity

58

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

Fig. 9. The vorticity contour for different We and at Re = 0. Contour levels are shown from −10 to 10 with increment of 1.

components decrease as expected. The maximum and minimum values of the velocity components at the corresponding centerlines are also given in Table 2 along with their locations. Here subscripts denote values pertaining to the minimum and maximum velocities. The results obtained at We = 0 compares well with those of a Newtonian fluid study [59]. Vertical position of umin shifts upward. Horizontal positions of the vmax and vmin , on the other hand, decrease as a function of We. In our calcu-

Fig. 10. Primary vortex center location as a function of We for Re = 0.

lations attainable limit value of We is 1 for creeping flow conditions. Table 2 also depicts that changes in horizontal and vertical velocities with respect to We become small after We ≥ 0.7. The cause of the little change on velocity field may be explained by momentum transfer from the lid to bulk. As elasticity increases momentum transfer may be confined to region close to the lid; hence effects of the lid do not go in to the deep of the cavity. Conversely as We decreases, momentum is transferred deeper in the cavity. The effect of We on the magnitude and location of the center of the primary eddy, values of the vorticity and first normal stress differences are listed in Table 3. The stream function value at the center of the primary eddy decreases as We increases. Other impact of the increasing We is the position of the primary eddy center moves slightly toward the upstream corner of the lid. Moreover, effect of the We is clearly depicted in Fig. 12, where upstream and downstream bottom corner streamlines with value of 0 are plotted at different We (the arrows point increasing We from 0 to the maximum value depicted in Table 3). Magnitudes and locations of the upstream secondary eddy (USE) and downstream secondary eddy (DSE) are also affected. Their height and width shrink with increasing We. Computed first normal stress difference, N1 , profiles are presented in Fig. 13 along the line y = 0.9989 as a function of x-direction at different We. Increasing elasticity leads to higher N1 due to the hoop stresses in the vicinity of the lid corner. In addition, at low We values, such as We = 0.1, N1 reaches its maximum at location of x = 0.04. After that it starts to relax to its minimum. Near the downstream corner, on the other hand, due to singularity effect its value starts to increase. At higher We N1 change with respect

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

59

Fig. 11. Profiles of computed (a) u along the line x = 0.5 and (b) v along the line y = 0.5 at Re = 0. Table 2 Horizontal minimum velocity, vertical minimum and maximum velocity through the centerline of the cavity at Re = 0. Data are obtained by M4. We

Reference

umin

ymin

vmax

xmax

vmin

xmin

0

Present [59]

−0.207738 −0.207754

0.535591 0.537600

0.184427 0.186273

0.210994 0.210500

−0.184427 −0.186273

0.789005 0.789400

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Present Present Present Present Present Present Present Present Present Present

−0.202928 −0.188620 −0.170218 −0.152492 −0.136846 −0.123454 −0.112087 −0.102455 −0.094280 −0.087298

0.539539 0.543485 0.551368 0.55530 0.559240 0.559240 0.559240 0.559240 0.555306 0.555306

0.180052 0.167173 0.151554 0.137284 0.125148 0.114997 0.106451 0.099216 0.093003 0.087691

0.210994 0.207631 0.207631 0.204281 0.204281 0.200946 0.197624 0.197624 0.194317 0.194317

−0.179086 −0.163109 −0.143441 −0.125829 −0.111347 −0.099653 −0.090138 −0.082329 −0.075843 −0.070379

0.789005 0.782238 0.775417 0.768543 0.761616 0.758113 0.751131 0.747611 0.744079 0.740534

Table 3 Properties of the primary vortex obtained by M4; the minimum stream function value, the vorticity value and first normal stress difference value and location of the center as a function of We for Re = 0. We

Reference

ω

N1

xmin

ymin

0

Present [59]

−0.100072 −0.100054

−3.212671 −3.220800

– –

0.500000 0.500000

0.765086 0.762600

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Present Present Present Present Present Present Present Present Present Present

−0.098264 −0.092494 −0.084793 −0.077339 −0.070754 −0.065109 −0.060281 −0.056149 −0.052595 −0.049501

−3.198255 −3.103465 −2.977111 −2.795894 −2.638557 −2.463854 −2.367557 −2.239667 −2.138173 −2.054843

0.524280 0.686163 0.660812 0.551703 0.477746 0.395511 0.381777 0.348610 0.330793 0.321571

0.500000 0.496042 0.488128 0.484171 0.480216 0.472309 0.468357 0.468357 0.464408 0.460460

0.768543 0.775411 0.785629 0.792368 0.799053 0.802375 0.808975 0.812254 0.815518 0.818768

min

Fig. 12. Comparison of the eddy sizes USE and DSE as a function of We in terms of contour level at Re = 0.

60

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

Table 4 Horizontal minimum velocity, vertical minimum and maximum velocity through the centerline of the cavity at Re = 100. Data are obtained by M4. We

Reference

umin

ymin

vmax

xmax

vmin

xmin

0

Present [59] [60] [61]

−0.213986 −0.213924 −0.210900 −0.214042

0.456514 0.459800 0.453100 0.458100

0.179536 0.180888 0.175270 0.179572

0.238383 0.235400 0.234400 0.237000

−0.253772 −0.256603 −0.245330 −0.253803

0.808975 0.812700 0.804700 0.810400

0.1 0.2 0.3 0.4 0.5 0.6 0.7

Present Present Present Present Present Present Present

−0.193175 −0.166454 −0.143382 −0.125055 −0.110449 −0.099034 −0.090093

0.488128 0.531642 0.578860 0.606138 0.613883 0.613883 0.606138

0.155133 0.125999 0.103111 0.090050 0.082311 0.076875 0.072529

0.234913 0.231456 0.221165 0.217761 0.214370 0.214370 0.210994

−0.211728 −0.156591 −0.115235 −0.093884 −0.082484 −0.075041 −0.069296

0.802375 0.785629 0.775417 0.771986 0.771986 0.771986 0.771986

Table 5 Properties of the primary vortex obtained by M4; the minimum stream function value and location of the center as a function of We for Re = 100.

Fig. 13. Profiles of computed normal stress differences as a function of We along the line of y = 0.9989 for Re = 0.

to x-direction exhibiting different behavior than that of lower We cases. The main distinction is that N1 does not relax to a minimum between lid corners at high We. The fluids with higher elasticity need possibly longer length than the current one. 4.2.4. Results at Re = 100 The impact of the non-linear inertial terms in the Cauchy stress equation gets amplified as Reynolds number increases and eventually flow becomes unstable when a critical Re is reached, i.e. Recrit ≈ 500 [46]. In this and following parts of the report, the effects of increasing Re (100 and 400) on the viscoelastic flow in the square

We

Reference

xmin

ymin

0

Present [59] [60] [62]

−0.103507 −0.103471 −0.103423 −0.102600

0.613883 0.618900 0.617200 0.617200

0.736977 0.740000 0.734400 0.734400

0.1 0.2 0.3 0.4 0.5 0.6 0.7

Present Present Present Present Present Present Present

−0.096136 −0.084883 −0.073872 −0.065387 −0.059021 −0.054099 −0.050140

0.621603 0.617747 0.606138 0.594478 0.578860 0.567099 0.559240

0.758134 0.785629 0.815518 0.831622 0.841106 0.847354 0.850455

min

lid are investigated. The velocity components through the vertical and horizontal centerlines of the cavity are compared with respect to the We in Fig. 14. Their values are tabulated in Table 4. As can be seen from both Fig. 14 and Table 4, as We increases the minimum value of the horizontal velocity component decreases in magnitude and its location moves closer to the lid. Also, extreme values of the vertical velocity components decrease and their minimum and maximum values approach to zero with increasing We. This trend depicts that as the elasticity increases recirculation flow inside a lid-driven cavity gets restricted. Moreover, effect of the increasing We on the hydrodynamic behavior of the flow is investigated with streamline and vorticity contours that are shown in Fig. 15 and Fig. 16, respectively. Streamline plots of Fig. 15 show that at high We primary vortex center moves in the upstream direction and streamlines near the downstream corners become more curved. Similar effects are observed for creeping flow. Locations of the primary vortex center as a function of We is presented in Fig. 17. Table 5 includes the intensities

Fig. 14. Profiles of computed (a) u along the line x = 0.5 and (b) v along the line y = 0.5 at Re = 100.

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

Fig. 15. The stream functions for different We at Re = 100. Contour levels shown for each plot are the same.

Fig. 16. The vorticity contour for different We and at Re = 100. Contour levels are shown from −10 to 10 with increment of 1.

61

62

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

Table 6 Horizontal minimum velocity, vertical minimum and maximum velocity through the centerlines of the cavity at Re = 400. Data are obtained by M4. We

Reference

umin

ymin

vmax

xmax

vmin

xmin

0

Present [59] [60]

−0.328491 −0.328375 −0.327260

0.280980 0.281500 0.281300

0.303630 0.304447 0.302030

0.224582 0.225300 0.226600

−0.453870 −0.456316 −0.449930

0.862706 0.862100 0.859400

0.1 0.2 0.3 0.4

Present Present Present Present

−0.296276 −0.250873 −0.144299 −0.078592

0.328927 0.448631 0.636934 0.682286

0.257662 0.160755 0.029944 0.035423

0.248868 0.299217 0.291890 0.2042817

−0.389682 −0.237622 −0.040831 −0.038701

0.837960 0.768543 0.744079 0.782238

Fig. 17. Primary vortex center location as a function of We at Re = 100.

Fig. 18. Profiles of computed normal stress differences as a function of We along the line of y = 0.9989 for Re = 100. Table 7 Properties of the primary vortex obtained by M4; the minimum stream function value and location of the center as a function of We for Re = 400.

of the primary eddy and location of the center as function of We. Their comparison with the literature is also given for Newtonian fluid. Results of the Re 100 indicate that attainable maximum Weissenberg number strongly depends on Reynolds number. At this Reynolds number maximum We values reached were 0.7. Other flow behaviors including DSE and USE sizes and N1 profile close to the lid exhibit similar trends to those in the Re = 0 case. For example in Fig. 18 the first normal stress difference, N1 , profiles along the line of y = 0.9989 up to We of 0.7 are shown. They match closely corresponding profiles at Re = 0 that are included in Fig. 13. The effect of convection seems to be not strong enough to impact the flow fields at Re = 100.

We

Reference

xmin

ymin

0

Present [59] [60]

−0.113935 −0.113897 −0.113909

0.555306 0.553600 0.554700

0.606138 0.607500 0.605500

0.1 0.2 0.3 0.4

Present Present Present Present

−0.103177 −0.085110 −0.060108 −0.045825

0.563172 0.586678 0.697103 0.751131

0.625454 0.686005 0.815518 0.886464

min

4.2.5. Results at Re = 400 In the final test, flow at Re = 400 is investigated. At this Reynolds number the maximum attainable We number at which a converged

Fig. 19. Profiles of computed (a) u along the line x = 0.5 and (b) v along the line y = 0.5 at Re = 400.

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

Fig. 20. The stream functions for different We at Re = 400. Contour levels shown for each plot are the same.

Fig. 21. The vorticity contour for different We and at Re = 400. Contour levels are shown from −10 to 10 with increment of 1.

63

64

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65

5. Conclusions In this study numerical calculations are carried out for steady laminar flow of Oldroyd-B fluid in a lid-driven square cavity geometry. The impacts of Reynolds and Weissenberg numbers are investigated by utilizing the SIMPLE algorithm on a non-uniform collocated grid system. Convective terms in the momentum equations are treated using central differences scheme and upwind approximation for the viscoelastic stresses. The following conclusions can then be drawn:

Fig. 22. Primary vortex center location as a function of We at Re = 400.

• More than one flow quantity should be considered to determine required mesh resolution for mesh independency in viscoelastic flows. • The convergence of the numerical scheme, in other words attainable maximum value of We is very sensitive to Reynolds number. • Upstream secondary eddy (USE) and downstream secondary eddy (DSE) sizes in lid-driven cavity are strongly affected by the value of We. For instance at creeping flow condition (Re = 0) and small Reynolds numbers (Re = 100) DSE and USE size decreases with increasing We. • At moderately large value of Reynolds number (Re = 400), DSE and USE merge and form second loop at the bottom of the cavity at We = 0.3. However when We is greater than 0.3, DSE and USE shrink and become similar to those at low Reynolds numbers. • Dimensionless magnitude of first normal stress differences, N1 , near the lid increases with increasing We for all Reynolds numbers. In this study several numerical experiments are performed using finite volume method based on collocated grids arrangement only. Therefore, we believe that systematic comparison of the staggered and collocated grid arrangements may be useful to clarify loss of convergence of viscoelastic fluid simulations in the future. Moreover we plan to approximate stress equations by high resolution schemes based on convective bounded criteria (CVC) such as MINMOD, OSHER, STOIC, COPLA, HOAB, VONOS in our future studies. Acknowledgement

Fig. 23. Profiles of computed normal stress differences as a function of We along to line of y = 0.9989 for Re = 400.

The authors acknowledge the financial support provided by The Scientific and Technological Research Council of Turkey through the research project 105M087. References

numerical solution was found is 0.4. Fig. 19 illustrates distribution of horizontal velocity along the vertical centerline and vertical velocity along the horizontal centerline. Their extreme values are listed in Table 6. From Fig. 19, as We increases the minimum value of the horizontal velocity component decreases in magnitude and its location moves closer to the lid. Vertical velocity, on the other hand is almost zero. Computed streamlines and vorticity contours are displayed in Figs. 20 and 21, respectively. Up to We = 0.3 DSE and USE both grow and at this We they merge and form a second loop at the cavity bottom. However, after We = 0.3 the second loop in the cavity disappears and DSE and USE begin to shrink with increasing. Fig. 22 shows location of the primary vortex center as a function of We. It is found that primary vortex shifts in downstream direction as opposed to the trend observed in Re = 0 and 100 simulations. Table 7 includes the intensities of the primary eddy with the various values of We. Their comparison with the literature is also given for Newtonian fluid. In Fig. 23, we demonstrate the graphs of normal stress differences along the lines of y = 0.9989.

[1] R.E. Evans, K. Walters, Further remarks on the lip-vortex mechanism of vortex enhancement in planar-contraction flows, J. Non-Newton. Fluid Mech. 32 (1989) 95–105. [2] R.E. Evans, K. Walters, Flow characteristic associated with abrupt changes in geometry in the case of highly elastic liquid, J. Non-Newton. Fluid Mech. 20 (1986) 11–29. [3] D.V. Boger, D.U. Hur, R.J. Binnington, Further observations of elastic effects in tubular entry flows, J. Non-Newton. Fluid Mech. 20 (1986) 31–49. [4] B. Pak, U.I. Cho, S.U.S. Choı, Separation and reattachment of non-Newtonian fluid flows in a sudden expansion pipe, J. Non-Newton. Fluid Mech. 37 (1990) 175–199. [5] Y.C. Zhou, B.S.V. Patnaik, D.C. Wan, G.W. Wei, DSC solution for flow in a staggered double lid driven cavity, Int. J. Numer. Meth. Eng. 57 (2003) 211–234. [6] M.S. Darwish, J.R. Whiteman, Numerical modeling of viscoelastic liquids using a finite-volume method, J. Non-Newton. Fluid Mech. 45 (1992) 311–337. [7] R. Keunings, On the high Weissenberg number problem, J. Non-Newton. Fluid Mech. 20 (1986) 209–226. [8] M. Fortin, A new approach for the FEM simulation of viscoelastic flows, J. NonNewton. Fluid Mech. 32 (1989) 295–310. [9] A.R. Davies, Numerical filtering and the high Weissenberg number problem, J. Non-Newton. Fluid Mech. 16 (1984) 195–209. [10] B. Debbaut, J.M. Marchal, M.J. Crochet, Numerical simulation of highly viscoelastic flows through an abrupt contraction, J. Non-Newton. Fluid Mech. 29 (1988) 119–146. [11] P. Saramito, Efficient simulation of nonlinear viscoelastic fluid flows, J. NonNewton. Fluid Mech. 60 (1995) 199–223.

K. Yapici et al. / J. Non-Newtonian Fluid Mech. 164 (2009) 51–65 [12] M.R. Apelian, R.C. Armstrong, R.A. Brown, Impact of the constitutive equation and singularity on the calculation of stick-slip flow: the modified upper-convected Maxwell model, J. Non-Newton. Fluid Mech. 27 (1988) 299–321. [13] G.P. Sasmal, A finite volume approach for calculation of viscoelastic flow through an abrupt axisymmetric contraction, J. Non-Newton. Fluid Mech. 56 (1995) 15–47. [14] D.D. Joseph, M. Renardy, J.-C. Saut, Hyperbolicity and change of type in the flow of viscoelastic fluids, Mech. Anal. 87 (1985) 213–251. [15] J.Y. Yoo, D.D. Joseph, Hyperbolicity and change of type in the flow of viscoelastic fluids through channels, J. Non-Newton. Fluid Mech. 19 (1985) 15–41. [16] T.B. Gasti, Steady flow of a non-Newtonian fluid through a contraction, J. Comput. Phys. 27 (1975) 42–70. [17] M.J. Crochet, G. Pilate, Plane flow of a fluid of second grade through a contraction, J. Non-Newton. Fluid Mech. 1 (1976) 247–258. [18] H.C. Choi, J.H. Song, J.Y. Yoo, Numerical simulation of planar contraction flow of a Gieskus fluid, J. Non-Newton. Fluid Mech. 29 (1988) 347–379. [19] P.W. Chang, T.W. Patten, B.A. Finlayson, Collocation and Galerkin finite element methods for viscoelastic fluid flow, Comput. Fluids 7 (1979) 267–283. [20] R.E. Gaidos, R. Darby, Numerical simulation and change in type in the developing flow of a nonlinear viscoelastic fluid, J. Non-Newton. Fluid Mech. 29 (1988) 59–79. [21] A.N. Beris, R.C. Armstrong, R.A. Brown, Spectral/finite element calculations of the flow of a Maxwell fluid between eccentric rotating cylinders, J. NonNewton. Fluid Mech. 22 (1987) 129–167. [22] J.Y. Yoo, Y. Na, A numerical study of the planar contraction flow of a viscoelastic fluid using the SIMPLER algorithm, J. Non-Newton. Fluid Mech. 39 (1991) 89–106. [23] S.-C. Xue, N. Phan-Thien, R.I. Tanner, Numerical study of secondary flows of viscoelastic fluid in straight pipes by an implicit finite volume method, J. NonNewton. Fluid Mech. 59 (1995) 191–213. [24] S.-C. Xue, N. Phan-Thien, R.I. Tanner, Three dimensional numerical simulations of viscoelastic flows through planar contractions, J. Non-Newton. Fluid Mech. 74 (1998) 195–245. [25] 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-Newton. Fluid Mech. 79 (1998) 1–43. [26] X. Huang, N. Phan-Thien, R.I. Tanner, Viscoelastic flow between eccentric rotating cylinders: unstructured control volume method, J. Non-Newton. Fluid Mech. 64 (1996) 71–92. [27] R.I. Tanner, S.-C. Xue, Computing transient flows with high elasticity, KoreaAust. Rheol. J. 14 (2002) 143–159. [28] H.A. Moatassime, D. Esselaoui, A finite volume approach for unsteady viscoelastic fluid flows, Int. J. Numer. Meth. Fluids 39 (2002) 939–959. [29] M. Peric, R. Kesser, G. Scheuerer, Comparison of finite volume numerical methods with staggered and collocated grids, Comput. Fluids 16 (1988) 389–403. [30] A.W. Date, Solution of Navier–Stokes equations on non-staggered grid, Int. J. Heat Mass Transfer 36 (1993) 1913–1922. [31] C.M. Rhie, W.L. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA J. 21 (11) (1983) 1525–1532. [32] S. Majumdar, Role of underrelaxation in momentum interpolation for calculation of flow with nonstaggered grids, Numer. Heat Transfer 13 (1988) 125–132. [33] B. Yu, W.-Q. Tau, J.-J. Wei, Y. Kawaguchi, T. Tagawa, H. Ozoe, Discussion on momentum interpolation method for calculated grids of incompressible flow, Numer. Heat Transfer B 42 (2002) 141–166. [34] M.A. Alves, F.T. Pinho, P.J. Oliveira, Effect of a high-resolution differencing scheme on finite-volume prediction of viscoelastic flows, J. Non-Newton. Fluid Mech. 93 (2000) 287–314. [35] M.A. Alves, F.T. Pinho, P.J. Oliveira, The flow of viscoelastic fluids past a cylinder: finite volume high-resolution methods, J. Non-Newton. Fluid Mech. 97 (2001) 207–232. [36] S.S. Edussuriya, A.J. Williams, C. Bailey, A cell-centred finite volume method for modeling viscoelastic flow, J. Non-Newton. Fluid Mech. 117 (2004) 47–61.

65

[37] D.B. Spalding, A novel finite difference formulation for differential expressions involving both first and second derivatives, Int. J. Numer. Methods Eng. 4 (1972) 551–559. [38] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York, 1980. [39] G.D. Raithby, Skew upwind differencing schemes for problems involving fluid flow, Comput. Methods Appl. Mech. Eng. 19 (1976) 153–164. [40] W. Shyy, S. Thakur, J. Wright, Second order upwind and central difference schemes for recirculating flow computation, AIAA J. 30 (1999) 923–932. [41] B.P. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comput. Methods Appl. Mech. Eng. 19 (1979) 59–98. [42] P.H. Gaskell, A.K.C. Lau, Curvature-compensated convective transport: SMART, a new boundedness-preserving transport algorithm, Int. J. Numer. Methods Fluids 8 (1988) 617–641. [43] T. Cochrane, K. Walters, M.F. Webster, On Newtonian and non-Newtonian flow in complex geometries, Phil. Trans. R. Soc. Lond. A 301 (163) (1981) 163–181. [44] M.A. Mendelson, P.-W. Yeh, R.A. Brown, R.C. Armstrong, Approximation error in finite element calculations of viscoelastic fluid flows, J. Non-Newton. Fluid Mech. 10 (1982) 31–54. [45] P. Pakdel, S.H. Spiegelberg, G.H. McKinley, Cavity flows of elastic liquids: twodimensional flows, Phys. Fluids 9 (11) (1997) 3123–3140. [46] P. Pakdel, G.H. McKinley, Cavity flows of elastic liquids: purely elastic instabilities, Phys. Fluids 10 (1998) 1058–1070. [47] A.M. Grillet, B. Yang, B. Khomami, E.S.G. Shaqfeh, Modelling of viscoelastic lid driven cavity flow using finite element simulations, J. Non-Newton. Fluid Mech. 88 (1999) 99–131. [48] A.M. Grillet, E.S.G. Shaqfeh, B. Khomami, Observations of elastic instabilities in lid-driven cavity flow, J. Non-Newton. Fluid Mech. 94 (2000) 15–35. [49] R. Fattal, R. Kupferman, Time dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation, J. NonNewton. Fluid Mech. 126 (2005) 23–37. [50] P.J. Oliveira, Viscoelastic lid-driven cavity flows, in: CMNE/CILAMCE 2007, Porto, Portugal (AMPTAC), 2007, pp. 1–16. [51] R.B. Bird, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, vol. 1, Fluid Mechanics, Wiley, USA, 1987. [52] Y. Na, J.Y. Yoo, A finite volume technique to simulate the flow of a viscoelastic fluid, Comput. Mech. 8 (1991) 43–55. [53] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Prentice Hall, 1995. [54] S.V. Patankar, D.B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Int. J. Heat Mass Transfer 15 (1972) 1787–1806. [55] J.P. Van Doormal, G.D. Raithby, Enhancement of the SIMPLE method for predicting incompressible flows, Numer. Heat Transfer 7 (1984) 147–163. [56] R.I. Issa, Solution of implicitly discretized fluid flow equation by operator—splitting, J. Comput. Phys. 62 (1986) 40–65. [57] M.F. Tomé, N. Mangiavacchi, J.A. Cuminato, A. Castelo, S. McKee, Afinite difference technique for simulating unsteady viscoelastic free surface flows, J. Non-Newton. Fluid Mech. 106 (2002) 61–106. [58] Q.-H. Deng, G.-F. Tang, Special treatment of pressure correction based on continuity conservation in a pressure based algorithm, Numer. Heat Transfer B 42 (2002) 73–92. [59] M. Sahin, R.G. Ovens, A novel fully implicit finite volume method applied to the lid-driven cavity problem. Part I. High Reynolds number flow calculations, Int. J. Numer. Fluids 42 (2003) 57–77. [60] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387–411. [61] O. Botella, R. Peyret, Benchmark spectral results on the lid-driven cavity flow, Comput. Fluids 27 (4) (1998) 421–433. [62] C.-H. Bruneau, C. Jouron, An efficient scheme for solving steady incompressible Navier–Stokes equations, J. Comput. Phys. 89 (2) (1990) 389–413.