- Email: [email protected]
A 3rd order upwind finite volume method for generalised Newtonian fluid flows
Advances in Engineering Software 36 (2005) 664–680 www.elsevier.com/locate/advengsoft
A 3rd order upwind finite volume method for generalised Newtonian fluid flows Panagiotis Neofytou* Environmental Research Lab., INT-RP, National Centre for Scientific Research ‘Demokritos’, Agia Paraskevi, 15310 Athens, Greece Received 12 October 2004; received in revised form 4 March 2005; accepted 4 March 2005 Available online 6 May 2005
Abstract The non-Newtonian effects in the flow of non-Newtonian fluids that can be modelled with generalised Newtonian constitutive equations are investigated using a numerical scheme based on the finite volume formulation. Collocated arrangement of variables is used and the pressure-correction method in conjunction with the SIMPLE scheme is applied. In order to avoid the diffusive effects of low order schemes, the approximation of the convection terms is carried out using the QUICK differencing scheme, which is of third order of accuracy. For modelling viscous non-Newtonian behaviour, the Power-Law, Quemada and the modified Bingham and Casson models are employed. Validation of the code is carried out upon comparison with numerical data available in the literature. Exploiting the lid-driven cavity flow a twofold investigation is carried out regarding the non-Newtonian effects first by using different models and secondly by using the shearthinning and shear-thickening attributes of the fluid. q 2005 Elsevier Ltd. All rights reserved. Keywords: Numerical simulation; Finite volume method; Generalised Newtonian models; QUICK
1. Introduction Many engineering problems require the investigation of both Newtonian and non-Newtonian flows. The categories of viscous and viscoelastic non-Newtonian fluids and the variety of related constitutive equations they provide cover applications such as petroleum, lubricants, food industry and blood flow. A computational fluid dynamics (CFD) approach to a non-Newtonian problem is generally more complicated than a Newtonian case due to the complexity introduced by the constitutive equation of the employed model in the diffusion terms of the Navier–Stokes equations. Different numerical methods have been employed in the past in order to simulate viscous and viscoplastic fluid flows, the most widely known simulations of which are carried out using the Power-Law and Bingham models, respectively. The finite element method (FEM) is used in several studies such as the one by Syrja¨la¨ [24], which investigates the fully developed laminar flow of a Power-Law non-Newtonian fluid in * Tel.: C30 210 6503 403; fax: C30 210 6525 004. E-mail address: [email protected]
0965-9978/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advengsoft.2005.03.011
a rectangular duct. The FEM was also utilised for simulating viscoplastic fluid flows in pipes [21] using the Bingham model, isothermal and non-isothermal fluid flow [2], liddriven cavity flow of viscoelastic fluids [8] and for investigation of the extrudate swell problem [4,10]. While the finite volume (FV) method is applied mostly in viscoelastic fluid flows [11,14,26] there are also recent studies where a FV approach has been used in order to study the instability of viscous and viscoplastic fluid flows in a channel with a sudden expansion [13] or to investigate the effects of viscous and viscoplastic pulsating flows in a channel with a stenosis [12]. Several undesirable numerical effects in FV methods such as artificial diffusion are induced by low order interpolation of the convection term of the Navier–Stokes equations. In order to overcome these effects, higher order interpolation schemes have been developed. Among these a very popular scheme is the QUICK (Quadratic Upwind Interpolation for Convective Kinematics) scheme [9], which is 3rd order accurate and therefore possesses higher order of accuracy than lower order schemes. This scheme has been applied successfully in terms of efficiency and accuracy to Newtonian fluid flow and heat transfer cases [19]. The effectiveness of the scheme has also been evaluated through comparison with the first and second order schemes in
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
the study by Xue et al. [26] where it is reported that with the QUICK scheme the artificial numerical diffusion is controlled and the use of the scheme is extended to solving viscoelastic fluid flows based on the Olroyd-B model. In the present study first a numerical method is presented where the QUICK scheme is applied in conjunction with the SIMPLE method for solving nonNewtonian viscous flows. For validating the method the lid-driven cavity flow problem is utilised and the results are compared with corresponding results available in the literature for cases of both Newtonian and non-Newtonian fluids. Secondly an investigation of the non-Newtonian effects in a lid-driven cavity flow is presented. The nonNewtonian models employed are viscous (Power-Law, Quemada) and viscoplastic (Bingham, Casson). The latter category is expressed following the Papanastasiou [15] approximation, which overcomes its discontinuous character by modifying it as a viscous-like model for all range of shear stress values. The study is carried out in terms of investigating the effects on the flow field of (i) nondimensional parameters in the aforementioned models and (ii) the shear-thickening or shear-thinning attributes of the fluid.
where i,j are either x or y. The viscosity coefficient h is constant for Newtonian flows, whereas it is a function of the shear rate g_ in the case of non-Newtonian flows. The shear rate can be written in tensor form as g_ Z 2D Z VV CVV T ;
where VZ(u,v), u and v are the velocity components in the x and y directions, respectively, and D is the deformation tensor. Hence, (3) can be also written as _ t Z hg;
(5)
where t is the shear-stress tensor. Following the principle of material objectivity, i.e. h should remain unchanged regardless of the frame of reference, the expression for g_ involves the second invariant of the rate of deformation tensor IID, that is pffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ Z 2trðD2 Þ: jgj (6) Therefore, in the case of non-Newtonian fluids the general relation between h and g_ is _ g: _ t Z hðjgjÞ
(7)
_ g_ t Z hðgÞ$ 2.1. Non-Newtonian flow equations The problem is governed by the two-dimensional incompressible flow equations. These are the continuity equation ð V$dS Z 0 (1) S
where VZ[u, v]T, and momentum equations ð ð ð ð v ru dU C ru V$dS ZK pix $dSC ðtxx ix Ctxy iy Þ$dS vt S
S
ð
ð
ð
ð
U
S
S
S
v rv dU C rv V$dS ZK piy $dS C ðtyx ix Ctyy iy Þ$dS vt (2b) for the x and y direction, respectively. dS equals to n$dS (n is the unit vector normal to the surface dS) and ix, iy are the unit vectors in the x and y directions, respectively. For modelling viscous non-Newtonian flows, extra attention should be given to the diffusion terms (DT), i.e. second terms of the RHS of Eq. (2a) and (2b). The shear stress t can be expressed in terms of the shear rate g_ as (3)
(8)
and then expressed in tensorial form, as in (7), so that _ which is referred to as ‘effective viscosity’, can be hðjgjÞ, specified. The effective viscosity is then used for the calculation of the diffusion terms in the Navier–Stokes equations. Substituting (5) into (2a) and (2b), and non-dimensionalising all the variables in the governing equations, the diffusion terms then become ð 1 DTx0 Z hðg_ xx ix C g_ xy iy Þ$dS (9) rUNl S
S
(2a)
tij Zhg_ ij ;
(4)
For two- or three-dimensional numerical simulations, any non-Newtonian constitutive equation in simple shear form should first be written as
2. Governing equations
U
665
DTy0 Z
1 rUNl
ð
hðg_ yx ix C g_ yy iy Þ$dS
(10)
S
where the asterisk superscript denotes dimensionless quantities and r, UN and l are the reference density, velocity and length, respectively, utilised in the nondimensionalisation. The effective viscosity should be _ that is, expressed in terms of the non-dimensional jgj, _ Þ. According to the non-dimensionalisation procedure hðjgj _ Z jgj
_ jgj : UN=l
(11)
_ Þ can be decomposed as a product of Generally hðjgj _ , say f ðjgj _ Þ, constant part, say Co, and a function of jgj
666
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
Following (13), Re 0 , which in this case is referred to as RePL, is given by
so that _ Þ: _ Þ Z Co f ðjgj hðjgj
(12) RePL Z
Inserting (12) into (9) and (10), one can see that a term similar to a Reynolds number results in front of the diffusion terms; this term is defined by Re 0 Z
rUNl : C0
(13)
In the case where h is constant (i.e. Newtonian flow), _ ÞZ 1. C0Zh and f ðjgj 2.2. Non-Newtonian models Four non-Newtonian constitutive relationships have been employed for modelling the diffusive terms. These are: the Power-Law model, which is commonly used for modelling shear-thinning or shear-thickening behaviour of non-Newtonian fluids, the Bingham and Casson models, which are used for modelling viscoplastic fluids and the Quemada model which is used for modelling fluid suspensions. The Bingham and Casson models are expressed according to the Papanastasiou [15] approximation, due to which they obtain a purely viscous character. 2.2.1. Power-Law model The shear stress for a Power-Law model is given by t Z kg_ n :
(14)
where k is a constant relative to the properties of the fluid and n is the exponent according to which shearthinning or shear thickening behaviour is regarded when n!1 or nO1, respectively. Bringing (14) into the form of (8), yields t Z kg_ nK1 g_
(15)
and its tensorial form is written as _ nK1 g_ t Z kjgj
(16)
Comparing (16) with (7), one concludes that the effective viscosity takes the form _ Z kjgj _ nK1 : hðjgjÞ
_ Þ Z k hðjgj
nK1 UN _ nK1 : jgj nK1 l
(18)
According to (12), in the case of a Power-Law-model based flow one obtains C0 Z k
nK1 UN ; nK1 l
_ Þ Z jgj _ nK1 : f ðjgj
(19)
(20)
The characteristic parameter for a Power-Law-model based flow is RePL. 2.2.2. Bingham model Bingham constitutive equation is given by ( _ jtjO ty ; t Z ty C hNg; g_ Z 0;
jtj! ty ;
(21)
where ty is the yield stress and hN is the asymptotic _ viscosity (i.e. viscosity when g/N). Eq. (21) can also be written in the form of (8) as: ( _ jtjO ty ; t Z ty g_ C hN g; (22) g_ Z 0; jtj! ty : Therefore, its tensorial form is given by ty _ jtjO ty C hN g; tZ _ jgj : g_ Z 0; jtj! ty : 8 <
(23)
The difficulty in applying Bingham model in numerical schemes lies in its discontinuous character. Papanastasiou [15] proposed an alternative expression that overcomes this obstacle, that is, one equation for the whole range of shearstress values. According to this approach, (23) is re-written as
ty _ Kmjgj t Z hN C ð1 K e Þ g_ (24) _ jgj Comparing (24) with (7), one concludes that the effective viscosity is ty _ Z hN C (25) hðjgjÞ 1 K eKmjgj_ : _ jgj The dimensionless expression of (25) according to (11) is
0 Bn _ Þ Z hN 1 C 1 K eKm jgj_ hðjgj (26) _ jgj where m 0 Z ðm=ðl=UNÞÞ and Bn is the Bingham number defined by
(17)
Using (11), the dimensionless form of (17) is given by
rln : nK2 kUN
Bn Z
ty l : hNUN
(27)
A critical value of the exponent m 0 in (27) exists, that depends on the magnitude of the Bingham number beyond which results are practically unaffected by m 0 [22]. According to (12), in the case of a Bingham-model-based flow 0 Bn _ Þ Z 1 C 1 K eKm jgj_ ; C0 Z hN; (28) f ðjgj _ jgj
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
and taking into account (13), Re 0 , which in this case is referred to as ReBI, is defined as ReBI Z
rUNl hN
(29)
The characteristic parameters for a Bingham-model-based flow are ReBI and Bn. 2.2.3. Casson model The constitutive equation for a Casson model is ( pffiffiffi pffiffiffiffiffi pffiffiffiffiffiffiffiffiffi t Z ty C hNg_ ; jtjO ty ; g_ Z 0;
jtj! ty ;
(30)
where ty is the yield stress and hN is the asymptotic viscosity. The form of (8) yields for (30): 8 pffiffiffiffiffiffi2 < t Z pffiffiffiffiffiffiffi _ jtjO ty ; ty g_ C hN g; (31) : g_ Z 0; jtj! t : y
8 <
Therefore, its tensorial form is qffiffiffiffiffi pffiffiffiffiffiffi2 ty _ tZ _ C hN g; jtjO ty jgj
: g_ Z 0;
jtj! ty :
(32)
The Papanastasiou approximation is also used to obtain a continuous function for the Casson model yielding rffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 2 ty pffiffiffiffiffiffi _ K mjgj t Z hN C g_ (33) 1 Ke _ jgj This equation has been found [17] to approach Casson’s equation very satisfactorily for mO100. Comparing (33) with (7), one concludes that the effective viscosity is rffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 2 ty pffiffiffiffiffiffi _ K mjgj : (34) hðjgjÞ Z hN C 1 Ke _ jgj
667
2.2.4. Quemada model The Quemada model [18] was originally proposed to describe the behaviour of concentrated disperse systems. The shear stress is given by pffiffiffiffiffiffiffiffiffi K2 _ g_ 1 k0 C kN g= pffiffiffiffiffiffiffiffiffi c 4 _ g; (37) t Z hF 1 K 2 1 C g= _ g_ c where hF is the viscosity of the suspending medium and f is the volumetric concentration of the suspension. Parameters g_ c , kN and k0 are to account for additional viscous attributes of the suspension. Pertinent to previous studies [12,13], the values used for these parameters are fZ0.45, kNZ2.07 and koZ4.33. Eq. (37) is already written in the form of (8). Its tensorial form would be pffiffiffiffiffiffiffiffiffiffiffiffi !K2 _ g_ 1 k0 C kN jgj= pffiffiffiffiffiffiffiffiffiffiffiffi c 4 t Z hF 1 K g_ (38) 2 1 C jgj= _ g_ c Therefore the effective viscosity according to (7) is pffiffiffiffiffiffiffiffiffiffiffiffi !K2 _ g_ 1 k0 C kN jgj= pffiffiffiffiffiffiffiffiffiffiffiffi c 4 _ Z hF 1 K : (39) hðjgjÞ 2 1 C jgj= _ g_ c The dimensionless expression of (39) according to (11) is !K2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ =g_ c 1 k0 C kN jgj pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ Þ Z hF 1 K ; (40) hðjgj 4 2 1 C jgj _ =g_ c where g_ c Z
g_ c : UN=l
(41)
According to (12), in the case of a Quemada-modelbased flow, !K2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ =g_ c 1 k0 C kN jgj pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ Þ Z 1K C0 Z hF ; f ðjgj 4 2 1 C jgj _ =g_ c (42) 0
The dimensionless expression of (34) according to (11) is sffiffiffiffiffiffiffiffiffi " # pffiffiffiffiffiffiffiffiffiffi 2 Bn _ K m 0 jgj _ Þ Z hN 1 C hðjgj (35) 1 Ke _ jgj where m 0 Z ðm=ðl=UNÞÞ and Bn is the Bingham number defined in (27). According to (12), in the case of a Cassonmodel-based flow sffiffiffiffiffiffiffiffiffi " # pffiffiffiffiffiffiffiffiffiffi 2 Bn _ K m 0 jgj _ Þ Z 1C C0 Z hN; f ðjgj ; 1 Ke _ jgj (36) and taking into account (13), Re 0 , which in this case is referred to as ReCA, is defined equally as ReBI in (29). The characteristic parameters for a Casson-model-based flow are ReCA and Bn.
Taking into account (13), Re , which in this case is referred to as ReQU, is given by ReQU Z
rUNl : hF
(43)
The characteristic parameters for a Quemada-model based-flow are ReQU and g_ c .
3. Numerical model The first step in applying the finite volume method is the grid generation, which is to divide the computational domain into control volumes. Each control volume is a 4-faced or a 6- faced cell for a 2D domain or a 3D domain, respectively. Then the governing equations are written in discretised form or in other words expressed as a function of
668
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
their variables on discrete points (nodes) of the domain. For a collocated arrangement of variables the nodes are located at the cell centres. The derivatives employed in the equations are approximated by central or upwind expressions. Finally the discretised equations are required to satisfy the momentum and mass balance at every control volume via an iterative solution algorithm, which in the present case is the Semi Implicit Method for PressureLinked Equations (SIMPLE) [16]. 3.1. Discretisation A typical computational molecule for a collocated arrangement of variables is shown in Fig. 1 with a central Control Volume (CV) with centre P and its neighbouring CVs with centres E, N, S and W. The four cell faces of the central CV are notated accordingly as e, n, s and w. For the discretisation procedure, the surface or volume integrals in Eqs. (1),(2a) and (2b) are applied at the perimeter or volume of every CV and furthermore the surface integrals may be split into four CV face integrals. The analysis is focused on the expression of Eqs. (1) and (2a) on face ‘e’ (Fig. 1) unless stated otherwise whereas the corresponding expressions on the other faces for these equations as well as (2b) are assumed to be treated in the same way. Since the scheme is implicit, the equations require iterative solutions. Therefore the system of equations to be constructed is linear with unknowns the values of the velocities at the new iteration m. On the mth iteration, all non-linear terms are approximated by a product of a value from the preceding iteration (denoted with mK1) and a ‘new’ value (denoted with m). For the unknowns the superscript m will be dropped and the coefficients A as well as the known source terms Q of the system of equations are determined following a rationale similar with other studies [6,14].
3.1.1. Mass conservation For the conservation properties dictated by Eq. (1) the discretisation procedure yields ð V$dS Z m_ e zuemK1 Se (44) Se
where umK1 can be calculated by linear interpolation e between umK1 and umK1 P E . The discretisation of the terms of the momentum equation is carried out based on its non-dimensional form, which for (2a) is ð ð ð v udU C uV$dS Z K pix $dS vt S S U ð 1 vu vv C h Vu C ix C iy $dS ð45Þ rUNl vx vx S
3.1.2. Unsteady term The first term of Eq. (45) takes the following discretised form. ð v DU ðu K unP Þ u dU z (46) vt Dt P U
where n denotes the previous time step and the corresponding value is treated as source. 3.1.3. Convection term For the convection term of Eq. (45) an interpolation scheme should be used for the expression of u on the cell face. Here the QUICK [9] scheme is used which is third order accurate and therefore reduces the problem of numerical diffusion compared to the Upwind Differencing Scheme (UDS), the expression of which for the cell face ‘e’ is ( uP if m_ e O 0 UDS ue Z : (47) uE if m_ e ! 0 Furthermore it is more accurate than the Central Differencing Scheme (CDS), the expression of which for the cell face ‘e’ is
N
uCDS Z uE le C uP ð1 K le Þ e
w
where x K xP : le Z e xE K xP
ne n
W
w
P
e
E
s sw
(48)
se S
Fig. 1. Topology of a grid with collocated arrangement of variables.
The QUICK scheme employs three nodal values of the interpolated variable and therefore an extended computational molecule is assumed (Fig. 2). For the cell face ‘e’ the expression yields uQUICK e ( uP CC1CðuP KuW ÞmK1 CC2CðuE KuP ÞmK1 if m_ e O0 Z uP CC1KðuE KuEE ÞmK1 CC2KðuP KuE ÞmK1 if m_ e !0 (49)
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
669
v from the previous iteration at the four adjacent nodes. The terms in brackets (superscript mK1) are treated as source terms. In view of the aforementioned approximations, Eq. (45) applied for the CV with centre P can be written in discretised form as X x Ai ui Z QxP ; i Z E; W; N; S (53) AxP uP C Fig. 2. Computational molecule for the QUICK scheme.
i
where Ai are the coefficients of the unknowns and QP are the source terms. The discretised form of the y-momentum Eq. (2b) is obtained in the same way yielding X y AyP vP C Ai vi Z QyP ; i Z E; W; N; S: (54)
where Dx2e Dx2e ; C1K Z ; 4ðDxw CDxe ÞDxw 4ðDxe CDxee ÞDxee 2Dxw CDxe 2Dxee CDxe ; C2K Z : C2C Z 4ðDxw CDxe Þ 4ðDxe CDxee Þ
C1C Z
ð
i
Therefore the discretised form of the convection term is ru V$dSzruQUICK m_ e e
(50)
Se
where the first term (uP) in (49) is calculated implicitly whereas the remaining terms (superscript mK1) are treated as source terms. 3.1.4. Pressure term The pressure term of the x-momentum equation CV is discretised as ð K pix $dS zKðpe Se K pw Sw ÞmK1 (51) S
and pmK1 can be calculated by linear where pmK1 e w mK1 interpolation between pmK1 and pmK1 P E , and between pP mK1 and pW , respectively. This term is treated as source.
3.2. Solution method For solving the flow and pressure field the SIMPLE method is applied. In the first step of the method a pressure (pmK1) and a velocity field (umK1, vmK1) are guessed. Then the discretised expressions (53) and (54) of the xmomentum and y-momentum equations, respectively, are solved using the guessed field to yield the updated velocity field values u*, v*. These are nodal values and therefore the values on the cell faces are also needed. Since a collocated grid is employed, a special interpolation scheme [20] needs to be applied in order to avoid oscillating distribution of pressure. According to this scheme, the velocity on the cell face should be directly linked to the pressure values at the adjacent nodes, therefore for the cell face ‘e’ one obtains 2 3mK1 ð ð ue Z u e K Ce 4 pix $dS K pix $dS 5 (55) S
3.1.5. Diffusion term The discretised expression of the diffusion term of the xmomentum equation corresponding to the whole CV is ð vu vv vu ðtxx ix Ctxy iy Þ$dS Z 2 h Se C h Ch S vx e vx vy n n S vu vv vu u KuP K 2 h S C h Ch S zhe E S vx w w vx vy s s xE KxP e u KuP u KuW u KuS Chn N Sn Khw P Sw Khs P S x KxP xP KxW xP KxS s N u KuP v Kvnw u KuW C he E S Chn ne S Khw P S xE KxP e xne Kxnw n xP KxW w mK1 v Kvsw Khs se S ð52Þ xse Kxsw s The values of h are calculated from (12) using the values of u and v from the previous iteration. vne, vnw, vse, and vsw are calculated each by linear interpolation from the values of
S
Ðwhere ue is the interpolated value between uP and uE , S pix $dS is the interpolated value between the pressure terms involved in Ðthe x-momentum conservation at CVs P and E (Fig. 1) and S pix $dS is the value employing directly the pressure values at the nodes P and E (i.e. Ð S pix $dS zðpE K pP ÞSe ). The coefficient C accounts for geometric properties of the CVs P and E. The values u*, v* do not generally satisfy the continuity equation and therefore have to be corrected as follows
um Zu C u 0 ;
vm Zv C v 0 :
(56)
This is possible if the pressure terms are also corrected as pm Z pmK1 C p 0 . The expression of um e according to (55) is 2 3m ð ð m 4 Z u K C pi $dS K pix $dS 5 (57) um e x e e S
S
Subtracting (55) from (57) yields
670
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
ue0 Z u~ e0 K Ce ðpE0 K pP0 ÞSe
(58)
Moving lid u* =1
where e u~e0 Z u m e Ku 20
ð
1m
0
ð
1mK1 3
C Ce 4@ pix $dS A K @ pix $dS A S
5:
(59)
S
The corrected velocities are required to satisfy the continuity equation. If the mass flux through cell face ‘e’ is m_ e Z um e Se
(60)
then similar expression for the other cell faces of CV P results in X (61) m_ c Z Dm_ C Dm_ 0 Z 0; c Z e; w; n; s:
x y Fig. 3. Geometry of the cavity with the moving lid.
c
where X yc Sc Z Dm_ ; c
X
Similarly underrelaxation for u and v should also be used so that yc0 Sc Z Dm_ 0 ;
c Z e; w; n; s:
c
(62) where yZu when cZe, w and yZv when cZn, s. Substituting (62) in (61) with the terms yc0 as in (58) yields X X y~c0 Sc ; Ai pi0 Z KDm_ K AP pP0 C c i (63) i Z E; W; N; S; c Z e; w; n; s: where A are the coefficients of the unknown variables resulting from the discretised expressions for the pressure corrections. The last terms on the RHS of (63) are omitted because they lead to an extended computational molecule involving pressure correction terms of more nodes than only P, E, W, N and S and also because they involve the velocity corrections which are not yet known. That leads to the pressure-correction equation X AP pP0 C Ai pi0 Z KDm_ ; i Z E; W; N; S: (64) i
Applying (64) on all CVs of the domain results to a system of equations, the solutions of which gives the pressure corrections on all nodes. These are used further to correct the velocity values, which will now satisfy the continuity equations. However, they do not satisfy the momentum equations, so another iteration must be performed using the solutions from the previous time step as initial guesses with the iteration procedure continuing until reaching the desired accuracy. It should be noted that due to neglecting the last term of the RHS of (63) the SIMPLE algorithm may lead to divergence and therefore underrelaxation should be used so that pm Z pmK1 C apUF p 0 ;
0! apUF ! 1:
(65)
um Z u C auUF u 0 ;
0! auUF ! 1:
(66)
vm Z v C avUF v 0 ;
0! avUF ! 1:
(67)
The solution of the systems of Eqs. (53), (54) and (64) is carried out using Stone’s [23] Strongly Implicit Procedure (SIP). 3.3. Boundary conditions The lid-driven cavity setup consists of a two-dimensional, 1!1 (in dimensionless units) cavity with a moving lid, the velocity of which in conjunction with the shear properties of the fluid cause a recirculation of the flow in the cavity (Fig. 3). The computational domain consists of a x!y space of [0,1]![0,1] in dimensionless units. The wall boundary is at xZ0, xZ1cy2[0,1] and at yZ0cx2[0,1]. The lid is at yZ1cx2[0,1] and moves with a steady velocity ulid which is used as reference for non-dimensionalisation ðulid Z 1Þ. The boundary conditions are as follows: At both the lid and the wall boundaries, the pressure values are derived by extrapolation from the inner nodes, whereas the dimensionless velocity value at the wall boundary and lid boundary are set equal to 0 and 1, respectively, following the no-slip condition. A reference value of p*Z1 is set at an inner node of the grid and is held constant throughout the computations.
4. Results and discussion Prior to proceeding with the validation it is essential to ensure about the grid independence of the computations. The aim is to achieve a maximum tolerance of around 1% change of a certain variable from one grid level to the more refined one. For these purpose three different uniform grids
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
671
1
1.0
ap=0.01 ap=0.05 ap=0.10
0.8
ap=0.20
0.1
ap=0.30 ap=0.40
y
Residual
0.6
1E-2
0.4
1E-3 0.2
120x120 grid 200x200 grid 300x300 grid
0.0
1E-4 0
-0.4
-0.2
0.0
0.2
0.4
0.6
0.8
100 200 300 400 500 600 700 800 900 1000 1100
1.0
Iterations
u Fig. 4. u-velocity component along the vertical centreline of the cavity for three different grid-resolutions (Newtonian flow, ReZ1000).
were used namely a 120!120, a 200!200, and a 300!300 grid. The results are shown in Fig. 4 for the u-velocity component along the vertical centreline for the Newtonian case with ReZ1000. Selecting as a criterion the value of umin, it can be seen that grid independence is achieved with the 200!200 grid (Table 1) and therefore based on this grid refinement study this grid is used for all later computations.
Fig. 5. Convergence history for different values of the underrelaxation factor for pressure apUF for Power-Law fluid flow at RePLZ100, nZ1.5.
RESv Z
NALL X! 1
! X y ! y ! Ai vi K QyP ! ; !AP vP C i
i Z E; W; N; S
I
(69)
RESp Z
NALL X! 1
! X ! ! Ai pi0 C Dm_ ! ; !AP pP0 C i
I
i Z E; W; N; S (70)
4.1. Code attributes The optimisation of convergence rate is an important aspect of the code attributes. Towards this a study as regards to the underrelaxation factors (UF) for p, u and v in Eqs. (65)–(67) is carried out. Computations are performed for different values of apUF with au;v UF Z 0:8 and for different p values of au;v with a Z 0:2 for the flow of a Power-Law UF UF fluid at RePLZ100 and nZ1.5. The convergence criterion is that the largest of residuals defined from Eqs. (53), (54) and (64) as
RESu Z
NALL X! 1
! X ! x ! Axi ui K QxP ! ; !AP uP C i
I
i Z E; W; N; S (68)
Table 1 Values of umin from the u-velocity profile along the vertical centreline of the cavity for three different grid-resolutions (Newtonian flow, ReZ1000) Mesh
umin
Change (%)
300!300 200!200 120!120
K0.3793 K0.3752 K0.3664
1 2.3 –
where the sum 1-NALL implies all control volumes of the domain, is to drop below 10K4. The results for the convergence history for different pressure and velocity UF are shown in Figs. 5 and 6, respectively. It can be seen that the convergence rate at the beginning of the computations is similar for values of apUF between 0.05–0.4 but becomes oscillatory for apUF R 0:3 as the computation progresses and therefore convergence is not achieved for the latter values (Fig. 5). In addition the convergence rate for apUF Z 0:01 is slower compared to the other values. Furthermore, from the convergence history of different au;v UF values it can be seen that the convergence rate is low for small values of au;v UF whereas it optimises for au;v UF Z 0:90 (Fig. 6(i)). It is worthwhile mentioning that for values close to au;v UF Z 0:6 (Fig. 6(ii)) the convergence rate becomes oscillatory and therefore convergence is not achieved for these values whereas for values just above 0.9 there is total divergence (Fig. 6(i)). In view of the above apUF Z 0:2 and au;v UF Z 0:9 are used for all subsequent computations. The validation of the numerical method presented above is carried out for the case of the steady lid-driven cavity flow for which results are available in the literature for both Newtonian and non-Newtonian fluids. The popularity of the lid-driven cavity flow case in the CFD community has established it as a benchmark case for code validation
672
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
(i)
(i) 1E+1
1.0
au=0.4 au=0.7 au=0.8
1
0.8
au=0.9
au=0.6
0.1
0.6
au=0.5
y
Residual
au=0.91
1E-2
0.4 Present results (Re=1000) Ghia et al. (Re=1000)
1E-3
0.2
Botella and Peyret (Re=1000) Present results (Re=400) Ghia et al. (Re=400)
1E-4
0.0 0
1000
2000
3000
4000
5000
6000
7000
8000
–0.4
–0.2
0.0
0.2
Iterations (ii)
0.4
0.6
0.8
1.0
u (ii)
1
0.4
au=0.59 au=0.595 au=0.60
0.2
au=0.61
0.1
au=0.615 0.0
v
Residual
au=0.62 1E-2
–0.2 Present results (Re=1000) –0.4
1E-3
Ghia et al. (Re=1000) Botella and Peyret (Re=1000) Present results (Re=400)
–0.6
Ghia et al. (Re=400)
1E-4 0
500
1000
1500
2000
2500
3000
3500
0.0
0.2
0.4
0.6
0.8
1.0
x
Iterations Fig. 6. (i) Convergence history for different values of the underrelaxation u;v factor for velocities aUF for Power-Law fluid flow at RePLZ100, nZ1.5; (ii) focus on values near au;v UF Z 0:6:
Fig. 7. Comparison between present results and the results by previous researchers [3,7] for Newtonian fluid flow. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
among many researchers [3,5,7]. The results for Newtonian flow from the present study encompass the u-velocity and v-velocity profiles along the vertical and horizontal centreline of the cavity, respectively, for ReZ400 and ReZ1000 for which detailed results are available in the literature [3,7]. The u-velocity profile is shown in Fig. 7(i) and it can be seen that it is in very good agreement with the results by Ghia et al. [7] in the case of ReZ400 and with the results by Ghia et al. [7] and Botella and Peyret [3] in the case of ReZ1000. The same yields for the v-velocity results (Fig. 7(ii)). The validation of the code regarding non-Newtonian computations is carried out considering the Power-Law
model for the flow of which results from previous researchers are available [2]. In accordance with this, calculations have been performed for RePLZ100 and different values of the Power-Law exponent n. Fig. 8 shows the results for the u-velocity along the vertical centreline (Fig. 8(i)) and the v-velocity along the horizontal centreline (Fig. 8(ii)) for nZ0.5 and nZ1.5 which are in very good agreement with the results by Bell and Surana [2]. Furthermore it can be seen that for nZ0.5 the velocity values are substantially lower than for nZ1.5 due to the low effective viscosity that the latter case exhibits around the centre of the cavity which allows for higher velocity gradients to occur.
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
673
Table 2 Comparison between the UDS, CDS and QUICK schemes: peak values for the u- and v- velocities at the vertical and horizontal centrelines, respectively, of a Power-Law fluid flow at RePLZ100 and nZ1.5 for different grid resolutions
(i) 1.0
0.8
Mesh
Scheme
umin
vmax
vmin
32!32
UDS CDS QUICK
K0.2237 K0.2379 K0.2387
0.2197 0.2262 0.2267
K0.2601 K0.2708 K0.2707
48!48
UDS CDS QUICK
K0.2285 K0.2385 K0.2389
0.2225 0.2268 0.2270
K0.2629 K0.2690 K0.2690
80!80
UDS CDS QUICK
K0.2327 K0.2390 K0.2391
0.2245 0.2273 0.2274
K0.2644 K0.2684 K0.2684
120!120
UDS CDS QUICK
K0.2349 K0.2391 K0.2392
0.2255 0.2274 0.2275
K0.2656 K0.2682 K0.2682
200!200
UDS CDS QUICK
K0.2367 K0.2392 K0.2393
0.2267 0.2277 0.2278
K0.2667 K0.2683 K0.2683
y
0.6
0.4
Present results (n=1.5) 0.2
Bell and Surana (n=1.5) Present results (n=0.5) Bell and Surana (n=0.5)
0.0 –0.4
–0.2
0.0
0.2
0.4
0.6
0.8
1.0
u (ii) 0.4 Present results (n=1.5) Bell and Surana (n=1.5) Present results (n=0.5)
v
0.2
Bell and Surana (n=0.5)
0.0
–0.2
–0.4 0.0
0.2
0.4
0.6
0.8
1.0
x Fig. 8. Comparison between present results and the results by previous researchers [1] for Power-Law fluid flow at RePLZ100 and for different values of Power-Law exponent n. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
The effect of the interpolation schemes on the nonNewtonian flow computations is examined for different grid resolutions and for the flow of a Power-Law fluid at RePLZ100 and nZ1.5. Table 2 shows the values of umin at the vertical centreline and vmax and vmin at the horizontal centreline derived from computations using the UDS, the CDS and the QUICK interpolation schemes. One can see that the refinement of the grid progressively lowers the difference between the three schemes. The difference between the results for the CDS and QUICK schemes for low grid-resolutions tends to agreement for high
grid-resolutions. The difference between changing from UDS to CDS is far more profound than changing from CDS to QUICK. In particular the difference in the values of umin between UDS and CDS for a 36!36 grid is 0.0142 but only 0.0008 between CDS and QUICK. For a 200! 200 grid there is still some small difference between the UDS and CDS whereas the CDS and QUICK are practically in agreement. Furthermore, the QUICK scheme appears more accurate at low resolutions compared to the CDS and UDS schemes because the results with this scheme for low grid-resolutions are closest to the results for high grid-resolutions and for which the results by all schemes tend to agreement. Finally, the UDS scheme due to its low order of accuracy shows a higher sensitivity to grid refinement compared to the other schemes as expected. It should be noted that higher order schemes do not necessarily yield accurate numerical results but may in some cases result in unwanted numerical instability. This is evident in solutions of moving-discontinuity problems where these schemes exhibit oscillatory behaviour in the vicinity of discontinuities [1,25]. 4.2. Non-Newtonian effects 4.2.1. Similarity in non-Newtonian flows Similarity in non-Newtonian flows of a certain model is achieved when all parameters derived from the nondimensionalisation procedure are the same. Furthermore the values of these properties determine the viscous behaviour of the model. Following the non-dimensionalisation procedure described in Section 2.2 one can conclude that the parameters upon which similarity depends are: RePL and n for the Power-Law model, Re and Bn for the Bingham and Casson models and ReQU and g_ c for the Quemada model.
674
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
Fig. 9. Streamlines for Bingham fluid flow at ReBIZ100 for (i) BnZ0.01, (ii) BnZ0.1 and (iii) BnZ1.
Therefore the non-Newtonian effects of the aforementioned models can be studied by investigating the flow effects caused by changing the values of these parameters. 4.2.2. Effects of Non-Newtonian parameters In order to investigate the effects of non-Newtonian parameters in a lid-driven cavity flow, three additional models were used namely the Bingham, Casson and Quemada models. It is important at this point to stress out that the first two models belong to the category of viscoplastic fluids but take the form of viscous shearthinning models when the Papanastasiou approximation is used whereas the Quemada model belongs originally to the category of viscous non-Newtonian models. The Reynolds numbers considered are ReBIZ100, ReCAZ100, ReQUZ 258.2 and were chosen so that for very high shear-stress values, all models tend to a Newtonian behaviour of ReZ100. The investigation was focused on the effects of the other characteristic parameter of every model namely the Bn number for the Bingham and Casson models and g_ c number for the Quemada flow. Therefore, three cases are examined for every model concerning different values of the aforementioned parameters the difference between which is of an order of magnitude. The latter was chosen so as to
give a more comprehensive and evident view of the trend associated with these parameters, each of which from this point onwards shall be called second characteristic parameter (SCP) of the model related to it. The streamline patterns are shown in Fig. 9 for the Bingham model, Fig. 10 for the Casson model and Fig. 11 for the Quemada model. One can see that the centre of the primary vortex is moved towards the lid and the magnitude of the secondary vortices is decreased when the SCP is increased. The latter occurs due to the weakening of the vortices induced by the increasing of the viscous behaviour of the model. The profile of the u- and v-velocities along the vertical and horizontal centreline, respectively, is shown in Figs. 12–14 for the Bingham, Casson and Quemada flows, respectively. The u-velocity maximum is lower and it takes place closer to the lid for higher values of the SCP of each model and in addition the velocity at the area near the bottom of the cavity is lower. This occurs due to the fact that the strengthening of the viscous attributes of the flow by increasing of the SCP of each model restricts the convection rate of the momentum from the lid to the bottom of the cavity thus leading to low velocity values at that area. The peak values of the v-velocity along the horizontal centreline
Fig. 10. Streamlines for Casson fluid flow at ReCAZ100 for (i) BnZ0.01, (ii) BnZ0.1 and (iii) BnZ1.
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
675
Fig. 11. Streamlines for Quemada fluid flow at ReQUZ258.5 for (i) g_ c Z 0:01, (ii) g_ c Z 0:1 and (iii) g_ c Z 1:
(i) 1.0
0.8
0.8
0.6
0.6
y
y
(i) 1.0
0.4
0.4
0.2
0.2
Bn=0.01
Bn=0.01 Bn=0.1
Bn=0.1
Bn=1
Bn=1 0.0
0.0 –0.4
–0.2
0.0
0.2
0.4
0.6
0.8
–0.4
1.0
–0.2
0.0
0.2
u
0.4
0.6
0.8
1.0
u (ii) 0.2
(ii) 0.2 Bn=0.01
Bn=0.01
Bn=0.1
Bn=0.1
Bn=1
0.1
Bn=1
0.1
0.0
v
v
0.0
–0.1
–0.1
–0.2
–0.2
–0.3
–0.3 0.0
0.2
0.4
0.6
0.8
1.0
x Fig. 12. Velocity profiles for Bingham fluid flow at ReBIZ100 and different values of Bn number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
0.0
0.2
0.4
x
0.6
0.8
1.0
Fig. 13. Velocity profiles for Casson fluid flow at ReCAZ100 and different values of Bn number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
676
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680 0.82
(i) 1.0
Bingham Casson 0.80
0.6
0.78
Quemada
c
y
y
0.8
Bn,gc = 1 Bn,gc = 0.1 Bn,gc = 0.01
0.76
0.4
0.2
0.74
gc=0.01 gc=0.1 gc=1
0.0
0.72 –0.4
–0.2
0.0
0.2
0.4
0.6
0.8
1.0
(ii) 0.2 gc=0.01 gc=0.1
v
0.0
–0.1
–0.2
–0.3 0.2
0.4
0.58
0.60
0.62
0.64
0.6
0.8
Fig. 15. Position of centre of main vortex for each model flow at ReBIZ ReCAZ100, ReQUZ258.2 and for different SCP.
the main vortex for the Casson and Quemada models compared to the Bingham model. By increasing further the SCP of each model the vortex centre in the Casson and Quemada cases moves further towards both the lid and the horizontal centreline in a very similar fashion whereas for the Bingham model and for this extent of SCP increasing, the relocation of the vortex centre is towards the right wall. All the aforementioned observations are also evident when comparing between the velocity profiles of every model. It can be seen that for the Bingham model the change in the velocity profile when the SCP is changed from 0.01 to 0.1 is less significant whereas for the Casson and Quemada models the change is more marked and is also very similar. Further increasing of the SCP induces a more marked change in the velocity profiles for all model cases.
gc=1
0.1
0.0
0.56
x
u
1.0
x Fig. 14. Velocity profiles for Quemada fluid flow at ReQUZ258.2 and different values of g_ c ðhgc Þ number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
are also reduced with increasing of the viscous properties of the flow that obstruct high velocity gradients from occurring. The effects of the different models can also be observed from the position of the centre of the main vortex shown in Fig. 15. It can be seen that increasing the SCP of every model causes the primary vortex to move towards the lid. This can also be linked to the aforementioned phenomenon of the u-velocity being lower and taking place closer to the lid for higher values of the SCP. The change of the SCP from 0.01 to 0.1 causes a more marked displacement of
4.2.3. Effects of shear-thinning and shear-thickening properties The shear-thinning and shear-thickening properties of the flow are simulated by using the Power-Law model with different values of the exponent n (Eq. (14)) namely nO1 for shear-thickening- and n!1 for shear-thinning behaviour. The profile of the u- and v-velocities along the vertical and horizontal centreline, respectively, for different values of RePL are shown in Figs. 16–18 for Newtonian (nZ1), shearthinning and shear-thickening behaviour, respectively. It can be seen that for all cases the increasing of RePL causes higher gradients in addition with increasing of the peak values for the u- and v-velocities and can be attributed to the gradual reducing of the viscous effects that occurs with increasing RePL. The effect on the corresponding displacement of the primary vortex is shown in Fig. 20(i). It can be seen that the centre of the main vortex shifts towards the lid with
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
677
(i) 1.0
(i) 1.0
Re=50 Re=50
0.8
Re=100
0.8
Re=100
Re=200
Re=200
Re=300
Re=300 0.6
Re=500
0.6
y
y
Re=500
0.4
0.4
0.2
0.2
0.0
0.0 –0.4
–0.2
0.0
0.2
0.4
0.6
0.8
–0.4
1.0
–0.2
0.0
0.2
0.4
0.6
0.8
1.0
u
u (ii) 0.2
(ii) 0.4
Re=50
Re=50
Re=100
Re=100 0.2
Re=200
Re=200
Re=300
Re=300
Re=500
0.0
Re=500
v
v
0.0
–0.2 –0.2
–0.4
–0.4
–0.6 0.0
0.2
0.4
0.6
0.8
0.0
1.0
Fig. 16. Velocity profiles for Newtonian fluid flow at different values of Re number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
decreasing RePL in all cases whereas for shear thinning fluids the displacement is more marked. It can also be seen that from one point with decreasing RePL the vortex centre not further moves towards the upper right edge but changes direction towards the left edge of the cavity. The effect of the extent of the shear-thinning and shearthickening attributes in the u- and v-profiles along the vertical and horizontal centrelines respectively is shown for RePLZ100 and different values of n in Fig. 19. By decreasing n the u and v-velocity peaks decrease whereas the centre of the main vortex moves towards the upper right corner of the cavity as seen in Fig. 20(ii).
0.2
0.4
0.6
0.8
1.0
x
x
Fig. 17. Velocity profiles for Power-Law fluid flow for nZ0.5 and different values of RePL number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
The effect of altering RePL on the wall shear stress can be investigated from the peak values of the dimensionless wall shear stress defined as tw Z
1 g_ n RePL w
(71)
on the lower wall of the cavity. These values are shown in Table 3. It can be seen that in contrast to the vortex displacement the peak value of tw is more sensitive to the changing of RePL for shear-thickening rather than for shearthinning fluids. For RePLZ50 the maximum of tw for nZ1.5
678
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
(i) 1.0
(i) 1.0
Re=50
0.8
0.8
Re=100
n=3
Re=200
n=2
Re=300 0.6
n=1.5
0.6
Re=500
y
y
n=1 n=0.75
0.4
0.4
0.2
0.2
n=0.5
0.0
0.0 –0.4
–0.2
0.0
0.2
0.4
0.6
0.8
-0.4
1.0
-0.2
0.0
0.2
0.4
0.6
0.8
1.0
u
u (ii) 0.4
(ii) 0.4 Re=50
n=3
Re=100
n=2 n=1.5
Re=200
0.2
Re=300
n=1
0.2
n=0.75
Re=500
n=0.5
v
v
0.0 0.0
–0.2
-0.2 –0.4
–0.6
-0.4 0.0
0.2
0.4
0.6
0.8
1.0
x
0.0
0.2
0.4
0.6
0.8
1.0
x
Fig. 18. Velocity profiles for Power-Law fluid flow for nZ1.5 and different values of RePL number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
Fig. 19. Velocity profiles for Power-Law fluid flow at RePLZ100 and different values of n. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
is markedly greater than for nZ0.5 but it ends up smaller for RePLZ500. It can also be seen that for high values of RePL the higher values of the peak of tw occurs neither for high nor for low values of n but for values close to nZ1.
arrangement of variables. The convection terms are discretised using the QUICK scheme in order to avoid the artificial diffusion caused by interpolation schemes of low order of accuracy. Non-Newtonian flow computations were performed using four different non-Newtonian models namely the Power-Law, Bingham, Casson and Quemada models. The code is validated against both Newtonian and non-Newtonian lid-driven cavity flow cases. Results from the current study were in very good agreement with results by previous researchers for both Newtonian and Power-Law fluid flow cases. From the parametric study regarding the SCP of the Bingham, Casson and Quemada models
5. Summary A numerical code for the computation of viscous nonNewtonian flows has been presented. The code is based on the pressure-correction method and a collocated
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
679
Table 3 Maximum of the tw distribution at the bottom of the cavity for flow of Power-Law model with nZ1.5, 1 and 0.5 and different values of RePL
(i) 1.0
0.8
0.6
y
Re=50 Re=100 Re=200 Re=300 Re=500
0.4
0.2
n=0.5 n=1 n=1.5
0.0 0.0
0.2
0.4
0.6
0.8
1.0
RePL
nZ1.5
nZ1
nZ0.5
50 100 200 300 500
0.0201 0.0102 0.0053 0.0037 0.0025
0.0142 0.0076 0.0048 0.0042 0.0041
0.0090 0.0051 0.0039 0.0037 0.0035
fluids the increasing of Re number causes higher velocity gradients and subsequent moving of the centre of the main vortex towards the lid; (ii) the more shear-thinning the fluid is, the closer to the upper right edge the centre of the main vortex takes place; (iii) the peak value of the dimensionless wall shear stress ðtw Þ at the lower wall of the cavity is more sensitive to the changing of Re for shear-thickening fluids and can be higher or lower than the corresponding value for shear-thinning fluids depending on the value of Re.
x (ii) 1.0
References 0.8
0.6
y
n=3 n=2 n=1.5 n=1 n=0.75 n=0.5
0.4
0.2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
x Fig. 20. Position of centre of main vortex for flow of Power Law model: (i) with nZ1.5, 1 and 0.5 and different values of RePL; (ii) at RePLZ100 and different values of n.
the following conclusions are drawn: (i) the increasing of the SCP is followed by shortening of the distance between the centre of the main vortex and the lid whereas the profile for both the u- and v-velocities is flattened; (ii) the behaviour of the Bingham model is less sensitive to the changing of its SCP than that of the Casson and Quemada models; (iii) the Casson and Quemada models exhibit very similar trends as to the changing of their SCP. From the parametric study regarding the shear-thickening- and shearthinning effects on the flow the following conclusions are drawn: (i) for both shear-thinning and shear-thickening
[1] Anderson JD. Computational fluid dynamics. Singapore: McGrawHill; 1995. [2] Bell BC, Surana KS. p-Version least squares finite element formulation for two-dimensional, incompressible, non-Newtonian, isothermal and non-isothermal flow. Int J Num Meth Fluids 1994;18: 127–62. [3] Botella O, Peyret R. Benchmark spectral solutions on the lid-driven cavity flow. Comput Fluids 1998;27:421–33. [4] Brasseur E, Fyrillas MM, Georgiou GC, Crochet MJ. The timedependent extrudate-swell problem of an Olroyd-B fluid with slip along the wall. J Rheol 1998;42:549–66. [5] Dailey LD, Pletcher RH. Evaluation of multigrid acceleration for preconditioned time-accurate Navier–Stokes algorithms. Comput Fluids 1996;25:791–811. [6] Ferziger JH, Peric M. Computational methods for fluid dynamics. Berlin: Springer; 1999. [7] Ghia U, Ghia KN, Shin CT. High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method. J Comput Phys 1982;48:387–411. [8] Grillet AM, Yang B, Khomami B, Shaqfeh ESG. Modelling of viscoelastic lid driven cavity flow using finite element simulations. J Non-Newt Fluid Mech 1999;88:99–131. [9] Leonard BP. A stable and accurate convective modelling procedure based on quadradic upstream interpolation. Comput Meth Appl Mech Engng 1979;19:59–98. [10] Mitsoulis E. Three dimensional non-Newtonian computations of extrudate swell with the finite element method. Comput Meth Appl Mech Engng 1999;180:333–44. [11] Mompean G. On predicting abrupt contraction flows with differential and algebraic viscoelastic models. Comput Fluids 2002;31:935–56. [12] Neofytou P, Drikakis D. Effects of blood models on flows through a stenosis. Int J Num Meth Fluids 2003;43:597–635. [13] Neofytou P, Drikakis D. Non-Newtonian flow instability in a channel with a sudden expansion. J Non-Newt Fluid Mech 2003; 111:127–50. [14] Oliveira PJ, Pinho FT, Pinto GA. Numerical simulation of non-linear elastic flows with a general collocated finite volume method. J NonNewt Fluid Mech 1998;79:1–43.
680
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
[15] Papanastasiou TC. Flow of materials with yield. J Rheology 1987;31: 385–404. [16] Patankar SV, Spalding DB. A calculation procedure for heat, mass and momentum transfer in three dimensional parabolic flows. Int J Heat Mass Transfer 1972;15:1787–806. [17] Pham TV, Mitsoulis E. Entry and Exit Flows of Casson Fluids. Can J Chem Engng 1994;72:1080–4. [18] Quemada D. Rheology of Concentrated Disperse systems III. General features of the proposed non-newtonian model. Comparison with experimental data. Rheologica Acta 1978;17:643–53. [19] Rahman MM, Miettinen A, Siikonen T, Modified SIMPLE. formulation on a collocated grid with an assessment of the simplified QUICK scheme. Numer Heat Transf, Part B 1996;30: 291–314. [20] Rhie CM, Chow WL. A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation. AIAA J 1983;21: 1525–32.
[21] Saramito P, Roquet N. An adaptive finite element method for viscoplastic fluid flows in pipes. Comput Meth Appl Mech Engng 2001;190:5391–412. [22] Smyrnaios DN, Tsamopoulos JA. Squeeze flow of Bingham plastics. J Non-Newt Fluid Mech 2001;100:165–90. [23] Stone HL. Iterative solution of implicit approximations of multidimensional partial differential equations. SIAM J Num Analysis 1968;5:530–58. [24] Syrja¨la¨ S. Finite-element analysis of fully developed laminar flow of Power-Law non-Newtonian fluid in a rectangular duct. Int Comm Heat Mass Transf 1995;22:549–57. [25] Tannehill JC, Anderson DA, Pletcher RH. Computational fluid mechanics and heat transfer. Washington, DC: Taylor and Francis; 1997. [26] Xue SC, Phan-Tien N, Tanner RI. Upwind with deferred correction (UPDC): an effective implementation of higher order convection schemes for implicit finite volume methods. J Non-Newt Fluid Mech 2002;108:1–24.
A 3rd order upwind finite volume method for generalised Newtonian fluid flows Panagiotis Neofytou* Environmental Research Lab., INT-RP, National Centre for Scientific Research ‘Demokritos’, Agia Paraskevi, 15310 Athens, Greece Received 12 October 2004; received in revised form 4 March 2005; accepted 4 March 2005 Available online 6 May 2005
Abstract The non-Newtonian effects in the flow of non-Newtonian fluids that can be modelled with generalised Newtonian constitutive equations are investigated using a numerical scheme based on the finite volume formulation. Collocated arrangement of variables is used and the pressure-correction method in conjunction with the SIMPLE scheme is applied. In order to avoid the diffusive effects of low order schemes, the approximation of the convection terms is carried out using the QUICK differencing scheme, which is of third order of accuracy. For modelling viscous non-Newtonian behaviour, the Power-Law, Quemada and the modified Bingham and Casson models are employed. Validation of the code is carried out upon comparison with numerical data available in the literature. Exploiting the lid-driven cavity flow a twofold investigation is carried out regarding the non-Newtonian effects first by using different models and secondly by using the shearthinning and shear-thickening attributes of the fluid. q 2005 Elsevier Ltd. All rights reserved. Keywords: Numerical simulation; Finite volume method; Generalised Newtonian models; QUICK
1. Introduction Many engineering problems require the investigation of both Newtonian and non-Newtonian flows. The categories of viscous and viscoelastic non-Newtonian fluids and the variety of related constitutive equations they provide cover applications such as petroleum, lubricants, food industry and blood flow. A computational fluid dynamics (CFD) approach to a non-Newtonian problem is generally more complicated than a Newtonian case due to the complexity introduced by the constitutive equation of the employed model in the diffusion terms of the Navier–Stokes equations. Different numerical methods have been employed in the past in order to simulate viscous and viscoplastic fluid flows, the most widely known simulations of which are carried out using the Power-Law and Bingham models, respectively. The finite element method (FEM) is used in several studies such as the one by Syrja¨la¨ [24], which investigates the fully developed laminar flow of a Power-Law non-Newtonian fluid in * Tel.: C30 210 6503 403; fax: C30 210 6525 004. E-mail address: [email protected]
0965-9978/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advengsoft.2005.03.011
a rectangular duct. The FEM was also utilised for simulating viscoplastic fluid flows in pipes [21] using the Bingham model, isothermal and non-isothermal fluid flow [2], liddriven cavity flow of viscoelastic fluids [8] and for investigation of the extrudate swell problem [4,10]. While the finite volume (FV) method is applied mostly in viscoelastic fluid flows [11,14,26] there are also recent studies where a FV approach has been used in order to study the instability of viscous and viscoplastic fluid flows in a channel with a sudden expansion [13] or to investigate the effects of viscous and viscoplastic pulsating flows in a channel with a stenosis [12]. Several undesirable numerical effects in FV methods such as artificial diffusion are induced by low order interpolation of the convection term of the Navier–Stokes equations. In order to overcome these effects, higher order interpolation schemes have been developed. Among these a very popular scheme is the QUICK (Quadratic Upwind Interpolation for Convective Kinematics) scheme [9], which is 3rd order accurate and therefore possesses higher order of accuracy than lower order schemes. This scheme has been applied successfully in terms of efficiency and accuracy to Newtonian fluid flow and heat transfer cases [19]. The effectiveness of the scheme has also been evaluated through comparison with the first and second order schemes in
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
the study by Xue et al. [26] where it is reported that with the QUICK scheme the artificial numerical diffusion is controlled and the use of the scheme is extended to solving viscoelastic fluid flows based on the Olroyd-B model. In the present study first a numerical method is presented where the QUICK scheme is applied in conjunction with the SIMPLE method for solving nonNewtonian viscous flows. For validating the method the lid-driven cavity flow problem is utilised and the results are compared with corresponding results available in the literature for cases of both Newtonian and non-Newtonian fluids. Secondly an investigation of the non-Newtonian effects in a lid-driven cavity flow is presented. The nonNewtonian models employed are viscous (Power-Law, Quemada) and viscoplastic (Bingham, Casson). The latter category is expressed following the Papanastasiou [15] approximation, which overcomes its discontinuous character by modifying it as a viscous-like model for all range of shear stress values. The study is carried out in terms of investigating the effects on the flow field of (i) nondimensional parameters in the aforementioned models and (ii) the shear-thickening or shear-thinning attributes of the fluid.
where i,j are either x or y. The viscosity coefficient h is constant for Newtonian flows, whereas it is a function of the shear rate g_ in the case of non-Newtonian flows. The shear rate can be written in tensor form as g_ Z 2D Z VV CVV T ;
where VZ(u,v), u and v are the velocity components in the x and y directions, respectively, and D is the deformation tensor. Hence, (3) can be also written as _ t Z hg;
(5)
where t is the shear-stress tensor. Following the principle of material objectivity, i.e. h should remain unchanged regardless of the frame of reference, the expression for g_ involves the second invariant of the rate of deformation tensor IID, that is pffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ Z 2trðD2 Þ: jgj (6) Therefore, in the case of non-Newtonian fluids the general relation between h and g_ is _ g: _ t Z hðjgjÞ
(7)
_ g_ t Z hðgÞ$ 2.1. Non-Newtonian flow equations The problem is governed by the two-dimensional incompressible flow equations. These are the continuity equation ð V$dS Z 0 (1) S
where VZ[u, v]T, and momentum equations ð ð ð ð v ru dU C ru V$dS ZK pix $dSC ðtxx ix Ctxy iy Þ$dS vt S
S
ð
ð
ð
ð
U
S
S
S
v rv dU C rv V$dS ZK piy $dS C ðtyx ix Ctyy iy Þ$dS vt (2b) for the x and y direction, respectively. dS equals to n$dS (n is the unit vector normal to the surface dS) and ix, iy are the unit vectors in the x and y directions, respectively. For modelling viscous non-Newtonian flows, extra attention should be given to the diffusion terms (DT), i.e. second terms of the RHS of Eq. (2a) and (2b). The shear stress t can be expressed in terms of the shear rate g_ as (3)
(8)
and then expressed in tensorial form, as in (7), so that _ which is referred to as ‘effective viscosity’, can be hðjgjÞ, specified. The effective viscosity is then used for the calculation of the diffusion terms in the Navier–Stokes equations. Substituting (5) into (2a) and (2b), and non-dimensionalising all the variables in the governing equations, the diffusion terms then become ð 1 DTx0 Z hðg_ xx ix C g_ xy iy Þ$dS (9) rUNl S
S
(2a)
tij Zhg_ ij ;
(4)
For two- or three-dimensional numerical simulations, any non-Newtonian constitutive equation in simple shear form should first be written as
2. Governing equations
U
665
DTy0 Z
1 rUNl
ð
hðg_ yx ix C g_ yy iy Þ$dS
(10)
S
where the asterisk superscript denotes dimensionless quantities and r, UN and l are the reference density, velocity and length, respectively, utilised in the nondimensionalisation. The effective viscosity should be _ that is, expressed in terms of the non-dimensional jgj, _ Þ. According to the non-dimensionalisation procedure hðjgj _ Z jgj
_ jgj : UN=l
(11)
_ Þ can be decomposed as a product of Generally hðjgj _ , say f ðjgj _ Þ, constant part, say Co, and a function of jgj
666
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
Following (13), Re 0 , which in this case is referred to as RePL, is given by
so that _ Þ: _ Þ Z Co f ðjgj hðjgj
(12) RePL Z
Inserting (12) into (9) and (10), one can see that a term similar to a Reynolds number results in front of the diffusion terms; this term is defined by Re 0 Z
rUNl : C0
(13)
In the case where h is constant (i.e. Newtonian flow), _ ÞZ 1. C0Zh and f ðjgj 2.2. Non-Newtonian models Four non-Newtonian constitutive relationships have been employed for modelling the diffusive terms. These are: the Power-Law model, which is commonly used for modelling shear-thinning or shear-thickening behaviour of non-Newtonian fluids, the Bingham and Casson models, which are used for modelling viscoplastic fluids and the Quemada model which is used for modelling fluid suspensions. The Bingham and Casson models are expressed according to the Papanastasiou [15] approximation, due to which they obtain a purely viscous character. 2.2.1. Power-Law model The shear stress for a Power-Law model is given by t Z kg_ n :
(14)
where k is a constant relative to the properties of the fluid and n is the exponent according to which shearthinning or shear thickening behaviour is regarded when n!1 or nO1, respectively. Bringing (14) into the form of (8), yields t Z kg_ nK1 g_
(15)
and its tensorial form is written as _ nK1 g_ t Z kjgj
(16)
Comparing (16) with (7), one concludes that the effective viscosity takes the form _ Z kjgj _ nK1 : hðjgjÞ
_ Þ Z k hðjgj
nK1 UN _ nK1 : jgj nK1 l
(18)
According to (12), in the case of a Power-Law-model based flow one obtains C0 Z k
nK1 UN ; nK1 l
_ Þ Z jgj _ nK1 : f ðjgj
(19)
(20)
The characteristic parameter for a Power-Law-model based flow is RePL. 2.2.2. Bingham model Bingham constitutive equation is given by ( _ jtjO ty ; t Z ty C hNg; g_ Z 0;
jtj! ty ;
(21)
where ty is the yield stress and hN is the asymptotic _ viscosity (i.e. viscosity when g/N). Eq. (21) can also be written in the form of (8) as: ( _ jtjO ty ; t Z ty g_ C hN g; (22) g_ Z 0; jtj! ty : Therefore, its tensorial form is given by ty _ jtjO ty C hN g; tZ _ jgj : g_ Z 0; jtj! ty : 8 <
(23)
The difficulty in applying Bingham model in numerical schemes lies in its discontinuous character. Papanastasiou [15] proposed an alternative expression that overcomes this obstacle, that is, one equation for the whole range of shearstress values. According to this approach, (23) is re-written as
ty _ Kmjgj t Z hN C ð1 K e Þ g_ (24) _ jgj Comparing (24) with (7), one concludes that the effective viscosity is ty _ Z hN C (25) hðjgjÞ 1 K eKmjgj_ : _ jgj The dimensionless expression of (25) according to (11) is
0 Bn _ Þ Z hN 1 C 1 K eKm jgj_ hðjgj (26) _ jgj where m 0 Z ðm=ðl=UNÞÞ and Bn is the Bingham number defined by
(17)
Using (11), the dimensionless form of (17) is given by
rln : nK2 kUN
Bn Z
ty l : hNUN
(27)
A critical value of the exponent m 0 in (27) exists, that depends on the magnitude of the Bingham number beyond which results are practically unaffected by m 0 [22]. According to (12), in the case of a Bingham-model-based flow 0 Bn _ Þ Z 1 C 1 K eKm jgj_ ; C0 Z hN; (28) f ðjgj _ jgj
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
and taking into account (13), Re 0 , which in this case is referred to as ReBI, is defined as ReBI Z
rUNl hN
(29)
The characteristic parameters for a Bingham-model-based flow are ReBI and Bn. 2.2.3. Casson model The constitutive equation for a Casson model is ( pffiffiffi pffiffiffiffiffi pffiffiffiffiffiffiffiffiffi t Z ty C hNg_ ; jtjO ty ; g_ Z 0;
jtj! ty ;
(30)
where ty is the yield stress and hN is the asymptotic viscosity. The form of (8) yields for (30): 8 pffiffiffiffiffiffi2 < t Z pffiffiffiffiffiffiffi _ jtjO ty ; ty g_ C hN g; (31) : g_ Z 0; jtj! t : y
8 <
Therefore, its tensorial form is qffiffiffiffiffi pffiffiffiffiffiffi2 ty _ tZ _ C hN g; jtjO ty jgj
: g_ Z 0;
jtj! ty :
(32)
The Papanastasiou approximation is also used to obtain a continuous function for the Casson model yielding rffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 2 ty pffiffiffiffiffiffi _ K mjgj t Z hN C g_ (33) 1 Ke _ jgj This equation has been found [17] to approach Casson’s equation very satisfactorily for mO100. Comparing (33) with (7), one concludes that the effective viscosity is rffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 2 ty pffiffiffiffiffiffi _ K mjgj : (34) hðjgjÞ Z hN C 1 Ke _ jgj
667
2.2.4. Quemada model The Quemada model [18] was originally proposed to describe the behaviour of concentrated disperse systems. The shear stress is given by pffiffiffiffiffiffiffiffiffi K2 _ g_ 1 k0 C kN g= pffiffiffiffiffiffiffiffiffi c 4 _ g; (37) t Z hF 1 K 2 1 C g= _ g_ c where hF is the viscosity of the suspending medium and f is the volumetric concentration of the suspension. Parameters g_ c , kN and k0 are to account for additional viscous attributes of the suspension. Pertinent to previous studies [12,13], the values used for these parameters are fZ0.45, kNZ2.07 and koZ4.33. Eq. (37) is already written in the form of (8). Its tensorial form would be pffiffiffiffiffiffiffiffiffiffiffiffi !K2 _ g_ 1 k0 C kN jgj= pffiffiffiffiffiffiffiffiffiffiffiffi c 4 t Z hF 1 K g_ (38) 2 1 C jgj= _ g_ c Therefore the effective viscosity according to (7) is pffiffiffiffiffiffiffiffiffiffiffiffi !K2 _ g_ 1 k0 C kN jgj= pffiffiffiffiffiffiffiffiffiffiffiffi c 4 _ Z hF 1 K : (39) hðjgjÞ 2 1 C jgj= _ g_ c The dimensionless expression of (39) according to (11) is !K2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ =g_ c 1 k0 C kN jgj pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ Þ Z hF 1 K ; (40) hðjgj 4 2 1 C jgj _ =g_ c where g_ c Z
g_ c : UN=l
(41)
According to (12), in the case of a Quemada-modelbased flow, !K2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ =g_ c 1 k0 C kN jgj pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ Þ Z 1K C0 Z hF ; f ðjgj 4 2 1 C jgj _ =g_ c (42) 0
The dimensionless expression of (34) according to (11) is sffiffiffiffiffiffiffiffiffi " # pffiffiffiffiffiffiffiffiffiffi 2 Bn _ K m 0 jgj _ Þ Z hN 1 C hðjgj (35) 1 Ke _ jgj where m 0 Z ðm=ðl=UNÞÞ and Bn is the Bingham number defined in (27). According to (12), in the case of a Cassonmodel-based flow sffiffiffiffiffiffiffiffiffi " # pffiffiffiffiffiffiffiffiffiffi 2 Bn _ K m 0 jgj _ Þ Z 1C C0 Z hN; f ðjgj ; 1 Ke _ jgj (36) and taking into account (13), Re 0 , which in this case is referred to as ReCA, is defined equally as ReBI in (29). The characteristic parameters for a Casson-model-based flow are ReCA and Bn.
Taking into account (13), Re , which in this case is referred to as ReQU, is given by ReQU Z
rUNl : hF
(43)
The characteristic parameters for a Quemada-model based-flow are ReQU and g_ c .
3. Numerical model The first step in applying the finite volume method is the grid generation, which is to divide the computational domain into control volumes. Each control volume is a 4-faced or a 6- faced cell for a 2D domain or a 3D domain, respectively. Then the governing equations are written in discretised form or in other words expressed as a function of
668
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
their variables on discrete points (nodes) of the domain. For a collocated arrangement of variables the nodes are located at the cell centres. The derivatives employed in the equations are approximated by central or upwind expressions. Finally the discretised equations are required to satisfy the momentum and mass balance at every control volume via an iterative solution algorithm, which in the present case is the Semi Implicit Method for PressureLinked Equations (SIMPLE) [16]. 3.1. Discretisation A typical computational molecule for a collocated arrangement of variables is shown in Fig. 1 with a central Control Volume (CV) with centre P and its neighbouring CVs with centres E, N, S and W. The four cell faces of the central CV are notated accordingly as e, n, s and w. For the discretisation procedure, the surface or volume integrals in Eqs. (1),(2a) and (2b) are applied at the perimeter or volume of every CV and furthermore the surface integrals may be split into four CV face integrals. The analysis is focused on the expression of Eqs. (1) and (2a) on face ‘e’ (Fig. 1) unless stated otherwise whereas the corresponding expressions on the other faces for these equations as well as (2b) are assumed to be treated in the same way. Since the scheme is implicit, the equations require iterative solutions. Therefore the system of equations to be constructed is linear with unknowns the values of the velocities at the new iteration m. On the mth iteration, all non-linear terms are approximated by a product of a value from the preceding iteration (denoted with mK1) and a ‘new’ value (denoted with m). For the unknowns the superscript m will be dropped and the coefficients A as well as the known source terms Q of the system of equations are determined following a rationale similar with other studies [6,14].
3.1.1. Mass conservation For the conservation properties dictated by Eq. (1) the discretisation procedure yields ð V$dS Z m_ e zuemK1 Se (44) Se
where umK1 can be calculated by linear interpolation e between umK1 and umK1 P E . The discretisation of the terms of the momentum equation is carried out based on its non-dimensional form, which for (2a) is ð ð ð v udU C uV$dS Z K pix $dS vt S S U ð 1 vu vv C h Vu C ix C iy $dS ð45Þ rUNl vx vx S
3.1.2. Unsteady term The first term of Eq. (45) takes the following discretised form. ð v DU ðu K unP Þ u dU z (46) vt Dt P U
where n denotes the previous time step and the corresponding value is treated as source. 3.1.3. Convection term For the convection term of Eq. (45) an interpolation scheme should be used for the expression of u on the cell face. Here the QUICK [9] scheme is used which is third order accurate and therefore reduces the problem of numerical diffusion compared to the Upwind Differencing Scheme (UDS), the expression of which for the cell face ‘e’ is ( uP if m_ e O 0 UDS ue Z : (47) uE if m_ e ! 0 Furthermore it is more accurate than the Central Differencing Scheme (CDS), the expression of which for the cell face ‘e’ is
N
uCDS Z uE le C uP ð1 K le Þ e
w
where x K xP : le Z e xE K xP
ne n
W
w
P
e
E
s sw
(48)
se S
Fig. 1. Topology of a grid with collocated arrangement of variables.
The QUICK scheme employs three nodal values of the interpolated variable and therefore an extended computational molecule is assumed (Fig. 2). For the cell face ‘e’ the expression yields uQUICK e ( uP CC1CðuP KuW ÞmK1 CC2CðuE KuP ÞmK1 if m_ e O0 Z uP CC1KðuE KuEE ÞmK1 CC2KðuP KuE ÞmK1 if m_ e !0 (49)
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
669
v from the previous iteration at the four adjacent nodes. The terms in brackets (superscript mK1) are treated as source terms. In view of the aforementioned approximations, Eq. (45) applied for the CV with centre P can be written in discretised form as X x Ai ui Z QxP ; i Z E; W; N; S (53) AxP uP C Fig. 2. Computational molecule for the QUICK scheme.
i
where Ai are the coefficients of the unknowns and QP are the source terms. The discretised form of the y-momentum Eq. (2b) is obtained in the same way yielding X y AyP vP C Ai vi Z QyP ; i Z E; W; N; S: (54)
where Dx2e Dx2e ; C1K Z ; 4ðDxw CDxe ÞDxw 4ðDxe CDxee ÞDxee 2Dxw CDxe 2Dxee CDxe ; C2K Z : C2C Z 4ðDxw CDxe Þ 4ðDxe CDxee Þ
C1C Z
ð
i
Therefore the discretised form of the convection term is ru V$dSzruQUICK m_ e e
(50)
Se
where the first term (uP) in (49) is calculated implicitly whereas the remaining terms (superscript mK1) are treated as source terms. 3.1.4. Pressure term The pressure term of the x-momentum equation CV is discretised as ð K pix $dS zKðpe Se K pw Sw ÞmK1 (51) S
and pmK1 can be calculated by linear where pmK1 e w mK1 interpolation between pmK1 and pmK1 P E , and between pP mK1 and pW , respectively. This term is treated as source.
3.2. Solution method For solving the flow and pressure field the SIMPLE method is applied. In the first step of the method a pressure (pmK1) and a velocity field (umK1, vmK1) are guessed. Then the discretised expressions (53) and (54) of the xmomentum and y-momentum equations, respectively, are solved using the guessed field to yield the updated velocity field values u*, v*. These are nodal values and therefore the values on the cell faces are also needed. Since a collocated grid is employed, a special interpolation scheme [20] needs to be applied in order to avoid oscillating distribution of pressure. According to this scheme, the velocity on the cell face should be directly linked to the pressure values at the adjacent nodes, therefore for the cell face ‘e’ one obtains 2 3mK1 ð ð ue Z u e K Ce 4 pix $dS K pix $dS 5 (55) S
3.1.5. Diffusion term The discretised expression of the diffusion term of the xmomentum equation corresponding to the whole CV is ð vu vv vu ðtxx ix Ctxy iy Þ$dS Z 2 h Se C h Ch S vx e vx vy n n S vu vv vu u KuP K 2 h S C h Ch S zhe E S vx w w vx vy s s xE KxP e u KuP u KuW u KuS Chn N Sn Khw P Sw Khs P S x KxP xP KxW xP KxS s N u KuP v Kvnw u KuW C he E S Chn ne S Khw P S xE KxP e xne Kxnw n xP KxW w mK1 v Kvsw Khs se S ð52Þ xse Kxsw s The values of h are calculated from (12) using the values of u and v from the previous iteration. vne, vnw, vse, and vsw are calculated each by linear interpolation from the values of
S
Ðwhere ue is the interpolated value between uP and uE , S pix $dS is the interpolated value between the pressure terms involved in Ðthe x-momentum conservation at CVs P and E (Fig. 1) and S pix $dS is the value employing directly the pressure values at the nodes P and E (i.e. Ð S pix $dS zðpE K pP ÞSe ). The coefficient C accounts for geometric properties of the CVs P and E. The values u*, v* do not generally satisfy the continuity equation and therefore have to be corrected as follows
um Zu C u 0 ;
vm Zv C v 0 :
(56)
This is possible if the pressure terms are also corrected as pm Z pmK1 C p 0 . The expression of um e according to (55) is 2 3m ð ð m 4 Z u K C pi $dS K pix $dS 5 (57) um e x e e S
S
Subtracting (55) from (57) yields
670
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
ue0 Z u~ e0 K Ce ðpE0 K pP0 ÞSe
(58)
Moving lid u* =1
where e u~e0 Z u m e Ku 20
ð
1m
0
ð
1mK1 3
C Ce 4@ pix $dS A K @ pix $dS A S
5:
(59)
S
The corrected velocities are required to satisfy the continuity equation. If the mass flux through cell face ‘e’ is m_ e Z um e Se
(60)
then similar expression for the other cell faces of CV P results in X (61) m_ c Z Dm_ C Dm_ 0 Z 0; c Z e; w; n; s:
x y Fig. 3. Geometry of the cavity with the moving lid.
c
where X yc Sc Z Dm_ ; c
X
Similarly underrelaxation for u and v should also be used so that yc0 Sc Z Dm_ 0 ;
c Z e; w; n; s:
c
(62) where yZu when cZe, w and yZv when cZn, s. Substituting (62) in (61) with the terms yc0 as in (58) yields X X y~c0 Sc ; Ai pi0 Z KDm_ K AP pP0 C c i (63) i Z E; W; N; S; c Z e; w; n; s: where A are the coefficients of the unknown variables resulting from the discretised expressions for the pressure corrections. The last terms on the RHS of (63) are omitted because they lead to an extended computational molecule involving pressure correction terms of more nodes than only P, E, W, N and S and also because they involve the velocity corrections which are not yet known. That leads to the pressure-correction equation X AP pP0 C Ai pi0 Z KDm_ ; i Z E; W; N; S: (64) i
Applying (64) on all CVs of the domain results to a system of equations, the solutions of which gives the pressure corrections on all nodes. These are used further to correct the velocity values, which will now satisfy the continuity equations. However, they do not satisfy the momentum equations, so another iteration must be performed using the solutions from the previous time step as initial guesses with the iteration procedure continuing until reaching the desired accuracy. It should be noted that due to neglecting the last term of the RHS of (63) the SIMPLE algorithm may lead to divergence and therefore underrelaxation should be used so that pm Z pmK1 C apUF p 0 ;
0! apUF ! 1:
(65)
um Z u C auUF u 0 ;
0! auUF ! 1:
(66)
vm Z v C avUF v 0 ;
0! avUF ! 1:
(67)
The solution of the systems of Eqs. (53), (54) and (64) is carried out using Stone’s [23] Strongly Implicit Procedure (SIP). 3.3. Boundary conditions The lid-driven cavity setup consists of a two-dimensional, 1!1 (in dimensionless units) cavity with a moving lid, the velocity of which in conjunction with the shear properties of the fluid cause a recirculation of the flow in the cavity (Fig. 3). The computational domain consists of a x!y space of [0,1]![0,1] in dimensionless units. The wall boundary is at xZ0, xZ1cy2[0,1] and at yZ0cx2[0,1]. The lid is at yZ1cx2[0,1] and moves with a steady velocity ulid which is used as reference for non-dimensionalisation ðulid Z 1Þ. The boundary conditions are as follows: At both the lid and the wall boundaries, the pressure values are derived by extrapolation from the inner nodes, whereas the dimensionless velocity value at the wall boundary and lid boundary are set equal to 0 and 1, respectively, following the no-slip condition. A reference value of p*Z1 is set at an inner node of the grid and is held constant throughout the computations.
4. Results and discussion Prior to proceeding with the validation it is essential to ensure about the grid independence of the computations. The aim is to achieve a maximum tolerance of around 1% change of a certain variable from one grid level to the more refined one. For these purpose three different uniform grids
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
671
1
1.0
ap=0.01 ap=0.05 ap=0.10
0.8
ap=0.20
0.1
ap=0.30 ap=0.40
y
Residual
0.6
1E-2
0.4
1E-3 0.2
120x120 grid 200x200 grid 300x300 grid
0.0
1E-4 0
-0.4
-0.2
0.0
0.2
0.4
0.6
0.8
100 200 300 400 500 600 700 800 900 1000 1100
1.0
Iterations
u Fig. 4. u-velocity component along the vertical centreline of the cavity for three different grid-resolutions (Newtonian flow, ReZ1000).
were used namely a 120!120, a 200!200, and a 300!300 grid. The results are shown in Fig. 4 for the u-velocity component along the vertical centreline for the Newtonian case with ReZ1000. Selecting as a criterion the value of umin, it can be seen that grid independence is achieved with the 200!200 grid (Table 1) and therefore based on this grid refinement study this grid is used for all later computations.
Fig. 5. Convergence history for different values of the underrelaxation factor for pressure apUF for Power-Law fluid flow at RePLZ100, nZ1.5.
RESv Z
NALL X! 1
! X y ! y ! Ai vi K QyP ! ; !AP vP C i
i Z E; W; N; S
I
(69)
RESp Z
NALL X! 1
! X ! ! Ai pi0 C Dm_ ! ; !AP pP0 C i
I
i Z E; W; N; S (70)
4.1. Code attributes The optimisation of convergence rate is an important aspect of the code attributes. Towards this a study as regards to the underrelaxation factors (UF) for p, u and v in Eqs. (65)–(67) is carried out. Computations are performed for different values of apUF with au;v UF Z 0:8 and for different p values of au;v with a Z 0:2 for the flow of a Power-Law UF UF fluid at RePLZ100 and nZ1.5. The convergence criterion is that the largest of residuals defined from Eqs. (53), (54) and (64) as
RESu Z
NALL X! 1
! X ! x ! Axi ui K QxP ! ; !AP uP C i
I
i Z E; W; N; S (68)
Table 1 Values of umin from the u-velocity profile along the vertical centreline of the cavity for three different grid-resolutions (Newtonian flow, ReZ1000) Mesh
umin
Change (%)
300!300 200!200 120!120
K0.3793 K0.3752 K0.3664
1 2.3 –
where the sum 1-NALL implies all control volumes of the domain, is to drop below 10K4. The results for the convergence history for different pressure and velocity UF are shown in Figs. 5 and 6, respectively. It can be seen that the convergence rate at the beginning of the computations is similar for values of apUF between 0.05–0.4 but becomes oscillatory for apUF R 0:3 as the computation progresses and therefore convergence is not achieved for the latter values (Fig. 5). In addition the convergence rate for apUF Z 0:01 is slower compared to the other values. Furthermore, from the convergence history of different au;v UF values it can be seen that the convergence rate is low for small values of au;v UF whereas it optimises for au;v UF Z 0:90 (Fig. 6(i)). It is worthwhile mentioning that for values close to au;v UF Z 0:6 (Fig. 6(ii)) the convergence rate becomes oscillatory and therefore convergence is not achieved for these values whereas for values just above 0.9 there is total divergence (Fig. 6(i)). In view of the above apUF Z 0:2 and au;v UF Z 0:9 are used for all subsequent computations. The validation of the numerical method presented above is carried out for the case of the steady lid-driven cavity flow for which results are available in the literature for both Newtonian and non-Newtonian fluids. The popularity of the lid-driven cavity flow case in the CFD community has established it as a benchmark case for code validation
672
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
(i)
(i) 1E+1
1.0
au=0.4 au=0.7 au=0.8
1
0.8
au=0.9
au=0.6
0.1
0.6
au=0.5
y
Residual
au=0.91
1E-2
0.4 Present results (Re=1000) Ghia et al. (Re=1000)
1E-3
0.2
Botella and Peyret (Re=1000) Present results (Re=400) Ghia et al. (Re=400)
1E-4
0.0 0
1000
2000
3000
4000
5000
6000
7000
8000
–0.4
–0.2
0.0
0.2
Iterations (ii)
0.4
0.6
0.8
1.0
u (ii)
1
0.4
au=0.59 au=0.595 au=0.60
0.2
au=0.61
0.1
au=0.615 0.0
v
Residual
au=0.62 1E-2
–0.2 Present results (Re=1000) –0.4
1E-3
Ghia et al. (Re=1000) Botella and Peyret (Re=1000) Present results (Re=400)
–0.6
Ghia et al. (Re=400)
1E-4 0
500
1000
1500
2000
2500
3000
3500
0.0
0.2
0.4
0.6
0.8
1.0
x
Iterations Fig. 6. (i) Convergence history for different values of the underrelaxation u;v factor for velocities aUF for Power-Law fluid flow at RePLZ100, nZ1.5; (ii) focus on values near au;v UF Z 0:6:
Fig. 7. Comparison between present results and the results by previous researchers [3,7] for Newtonian fluid flow. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
among many researchers [3,5,7]. The results for Newtonian flow from the present study encompass the u-velocity and v-velocity profiles along the vertical and horizontal centreline of the cavity, respectively, for ReZ400 and ReZ1000 for which detailed results are available in the literature [3,7]. The u-velocity profile is shown in Fig. 7(i) and it can be seen that it is in very good agreement with the results by Ghia et al. [7] in the case of ReZ400 and with the results by Ghia et al. [7] and Botella and Peyret [3] in the case of ReZ1000. The same yields for the v-velocity results (Fig. 7(ii)). The validation of the code regarding non-Newtonian computations is carried out considering the Power-Law
model for the flow of which results from previous researchers are available [2]. In accordance with this, calculations have been performed for RePLZ100 and different values of the Power-Law exponent n. Fig. 8 shows the results for the u-velocity along the vertical centreline (Fig. 8(i)) and the v-velocity along the horizontal centreline (Fig. 8(ii)) for nZ0.5 and nZ1.5 which are in very good agreement with the results by Bell and Surana [2]. Furthermore it can be seen that for nZ0.5 the velocity values are substantially lower than for nZ1.5 due to the low effective viscosity that the latter case exhibits around the centre of the cavity which allows for higher velocity gradients to occur.
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
673
Table 2 Comparison between the UDS, CDS and QUICK schemes: peak values for the u- and v- velocities at the vertical and horizontal centrelines, respectively, of a Power-Law fluid flow at RePLZ100 and nZ1.5 for different grid resolutions
(i) 1.0
0.8
Mesh
Scheme
umin
vmax
vmin
32!32
UDS CDS QUICK
K0.2237 K0.2379 K0.2387
0.2197 0.2262 0.2267
K0.2601 K0.2708 K0.2707
48!48
UDS CDS QUICK
K0.2285 K0.2385 K0.2389
0.2225 0.2268 0.2270
K0.2629 K0.2690 K0.2690
80!80
UDS CDS QUICK
K0.2327 K0.2390 K0.2391
0.2245 0.2273 0.2274
K0.2644 K0.2684 K0.2684
120!120
UDS CDS QUICK
K0.2349 K0.2391 K0.2392
0.2255 0.2274 0.2275
K0.2656 K0.2682 K0.2682
200!200
UDS CDS QUICK
K0.2367 K0.2392 K0.2393
0.2267 0.2277 0.2278
K0.2667 K0.2683 K0.2683
y
0.6
0.4
Present results (n=1.5) 0.2
Bell and Surana (n=1.5) Present results (n=0.5) Bell and Surana (n=0.5)
0.0 –0.4
–0.2
0.0
0.2
0.4
0.6
0.8
1.0
u (ii) 0.4 Present results (n=1.5) Bell and Surana (n=1.5) Present results (n=0.5)
v
0.2
Bell and Surana (n=0.5)
0.0
–0.2
–0.4 0.0
0.2
0.4
0.6
0.8
1.0
x Fig. 8. Comparison between present results and the results by previous researchers [1] for Power-Law fluid flow at RePLZ100 and for different values of Power-Law exponent n. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
The effect of the interpolation schemes on the nonNewtonian flow computations is examined for different grid resolutions and for the flow of a Power-Law fluid at RePLZ100 and nZ1.5. Table 2 shows the values of umin at the vertical centreline and vmax and vmin at the horizontal centreline derived from computations using the UDS, the CDS and the QUICK interpolation schemes. One can see that the refinement of the grid progressively lowers the difference between the three schemes. The difference between the results for the CDS and QUICK schemes for low grid-resolutions tends to agreement for high
grid-resolutions. The difference between changing from UDS to CDS is far more profound than changing from CDS to QUICK. In particular the difference in the values of umin between UDS and CDS for a 36!36 grid is 0.0142 but only 0.0008 between CDS and QUICK. For a 200! 200 grid there is still some small difference between the UDS and CDS whereas the CDS and QUICK are practically in agreement. Furthermore, the QUICK scheme appears more accurate at low resolutions compared to the CDS and UDS schemes because the results with this scheme for low grid-resolutions are closest to the results for high grid-resolutions and for which the results by all schemes tend to agreement. Finally, the UDS scheme due to its low order of accuracy shows a higher sensitivity to grid refinement compared to the other schemes as expected. It should be noted that higher order schemes do not necessarily yield accurate numerical results but may in some cases result in unwanted numerical instability. This is evident in solutions of moving-discontinuity problems where these schemes exhibit oscillatory behaviour in the vicinity of discontinuities [1,25]. 4.2. Non-Newtonian effects 4.2.1. Similarity in non-Newtonian flows Similarity in non-Newtonian flows of a certain model is achieved when all parameters derived from the nondimensionalisation procedure are the same. Furthermore the values of these properties determine the viscous behaviour of the model. Following the non-dimensionalisation procedure described in Section 2.2 one can conclude that the parameters upon which similarity depends are: RePL and n for the Power-Law model, Re and Bn for the Bingham and Casson models and ReQU and g_ c for the Quemada model.
674
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
Fig. 9. Streamlines for Bingham fluid flow at ReBIZ100 for (i) BnZ0.01, (ii) BnZ0.1 and (iii) BnZ1.
Therefore the non-Newtonian effects of the aforementioned models can be studied by investigating the flow effects caused by changing the values of these parameters. 4.2.2. Effects of Non-Newtonian parameters In order to investigate the effects of non-Newtonian parameters in a lid-driven cavity flow, three additional models were used namely the Bingham, Casson and Quemada models. It is important at this point to stress out that the first two models belong to the category of viscoplastic fluids but take the form of viscous shearthinning models when the Papanastasiou approximation is used whereas the Quemada model belongs originally to the category of viscous non-Newtonian models. The Reynolds numbers considered are ReBIZ100, ReCAZ100, ReQUZ 258.2 and were chosen so that for very high shear-stress values, all models tend to a Newtonian behaviour of ReZ100. The investigation was focused on the effects of the other characteristic parameter of every model namely the Bn number for the Bingham and Casson models and g_ c number for the Quemada flow. Therefore, three cases are examined for every model concerning different values of the aforementioned parameters the difference between which is of an order of magnitude. The latter was chosen so as to
give a more comprehensive and evident view of the trend associated with these parameters, each of which from this point onwards shall be called second characteristic parameter (SCP) of the model related to it. The streamline patterns are shown in Fig. 9 for the Bingham model, Fig. 10 for the Casson model and Fig. 11 for the Quemada model. One can see that the centre of the primary vortex is moved towards the lid and the magnitude of the secondary vortices is decreased when the SCP is increased. The latter occurs due to the weakening of the vortices induced by the increasing of the viscous behaviour of the model. The profile of the u- and v-velocities along the vertical and horizontal centreline, respectively, is shown in Figs. 12–14 for the Bingham, Casson and Quemada flows, respectively. The u-velocity maximum is lower and it takes place closer to the lid for higher values of the SCP of each model and in addition the velocity at the area near the bottom of the cavity is lower. This occurs due to the fact that the strengthening of the viscous attributes of the flow by increasing of the SCP of each model restricts the convection rate of the momentum from the lid to the bottom of the cavity thus leading to low velocity values at that area. The peak values of the v-velocity along the horizontal centreline
Fig. 10. Streamlines for Casson fluid flow at ReCAZ100 for (i) BnZ0.01, (ii) BnZ0.1 and (iii) BnZ1.
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
675
Fig. 11. Streamlines for Quemada fluid flow at ReQUZ258.5 for (i) g_ c Z 0:01, (ii) g_ c Z 0:1 and (iii) g_ c Z 1:
(i) 1.0
0.8
0.8
0.6
0.6
y
y
(i) 1.0
0.4
0.4
0.2
0.2
Bn=0.01
Bn=0.01 Bn=0.1
Bn=0.1
Bn=1
Bn=1 0.0
0.0 –0.4
–0.2
0.0
0.2
0.4
0.6
0.8
–0.4
1.0
–0.2
0.0
0.2
u
0.4
0.6
0.8
1.0
u (ii) 0.2
(ii) 0.2 Bn=0.01
Bn=0.01
Bn=0.1
Bn=0.1
Bn=1
0.1
Bn=1
0.1
0.0
v
v
0.0
–0.1
–0.1
–0.2
–0.2
–0.3
–0.3 0.0
0.2
0.4
0.6
0.8
1.0
x Fig. 12. Velocity profiles for Bingham fluid flow at ReBIZ100 and different values of Bn number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
0.0
0.2
0.4
x
0.6
0.8
1.0
Fig. 13. Velocity profiles for Casson fluid flow at ReCAZ100 and different values of Bn number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
676
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680 0.82
(i) 1.0
Bingham Casson 0.80
0.6
0.78
Quemada
c
y
y
0.8
Bn,gc = 1 Bn,gc = 0.1 Bn,gc = 0.01
0.76
0.4
0.2
0.74
gc=0.01 gc=0.1 gc=1
0.0
0.72 –0.4
–0.2
0.0
0.2
0.4
0.6
0.8
1.0
(ii) 0.2 gc=0.01 gc=0.1
v
0.0
–0.1
–0.2
–0.3 0.2
0.4
0.58
0.60
0.62
0.64
0.6
0.8
Fig. 15. Position of centre of main vortex for each model flow at ReBIZ ReCAZ100, ReQUZ258.2 and for different SCP.
the main vortex for the Casson and Quemada models compared to the Bingham model. By increasing further the SCP of each model the vortex centre in the Casson and Quemada cases moves further towards both the lid and the horizontal centreline in a very similar fashion whereas for the Bingham model and for this extent of SCP increasing, the relocation of the vortex centre is towards the right wall. All the aforementioned observations are also evident when comparing between the velocity profiles of every model. It can be seen that for the Bingham model the change in the velocity profile when the SCP is changed from 0.01 to 0.1 is less significant whereas for the Casson and Quemada models the change is more marked and is also very similar. Further increasing of the SCP induces a more marked change in the velocity profiles for all model cases.
gc=1
0.1
0.0
0.56
x
u
1.0
x Fig. 14. Velocity profiles for Quemada fluid flow at ReQUZ258.2 and different values of g_ c ðhgc Þ number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
are also reduced with increasing of the viscous properties of the flow that obstruct high velocity gradients from occurring. The effects of the different models can also be observed from the position of the centre of the main vortex shown in Fig. 15. It can be seen that increasing the SCP of every model causes the primary vortex to move towards the lid. This can also be linked to the aforementioned phenomenon of the u-velocity being lower and taking place closer to the lid for higher values of the SCP. The change of the SCP from 0.01 to 0.1 causes a more marked displacement of
4.2.3. Effects of shear-thinning and shear-thickening properties The shear-thinning and shear-thickening properties of the flow are simulated by using the Power-Law model with different values of the exponent n (Eq. (14)) namely nO1 for shear-thickening- and n!1 for shear-thinning behaviour. The profile of the u- and v-velocities along the vertical and horizontal centreline, respectively, for different values of RePL are shown in Figs. 16–18 for Newtonian (nZ1), shearthinning and shear-thickening behaviour, respectively. It can be seen that for all cases the increasing of RePL causes higher gradients in addition with increasing of the peak values for the u- and v-velocities and can be attributed to the gradual reducing of the viscous effects that occurs with increasing RePL. The effect on the corresponding displacement of the primary vortex is shown in Fig. 20(i). It can be seen that the centre of the main vortex shifts towards the lid with
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
677
(i) 1.0
(i) 1.0
Re=50 Re=50
0.8
Re=100
0.8
Re=100
Re=200
Re=200
Re=300
Re=300 0.6
Re=500
0.6
y
y
Re=500
0.4
0.4
0.2
0.2
0.0
0.0 –0.4
–0.2
0.0
0.2
0.4
0.6
0.8
–0.4
1.0
–0.2
0.0
0.2
0.4
0.6
0.8
1.0
u
u (ii) 0.2
(ii) 0.4
Re=50
Re=50
Re=100
Re=100 0.2
Re=200
Re=200
Re=300
Re=300
Re=500
0.0
Re=500
v
v
0.0
–0.2 –0.2
–0.4
–0.4
–0.6 0.0
0.2
0.4
0.6
0.8
0.0
1.0
Fig. 16. Velocity profiles for Newtonian fluid flow at different values of Re number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
decreasing RePL in all cases whereas for shear thinning fluids the displacement is more marked. It can also be seen that from one point with decreasing RePL the vortex centre not further moves towards the upper right edge but changes direction towards the left edge of the cavity. The effect of the extent of the shear-thinning and shearthickening attributes in the u- and v-profiles along the vertical and horizontal centrelines respectively is shown for RePLZ100 and different values of n in Fig. 19. By decreasing n the u and v-velocity peaks decrease whereas the centre of the main vortex moves towards the upper right corner of the cavity as seen in Fig. 20(ii).
0.2
0.4
0.6
0.8
1.0
x
x
Fig. 17. Velocity profiles for Power-Law fluid flow for nZ0.5 and different values of RePL number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
The effect of altering RePL on the wall shear stress can be investigated from the peak values of the dimensionless wall shear stress defined as tw Z
1 g_ n RePL w
(71)
on the lower wall of the cavity. These values are shown in Table 3. It can be seen that in contrast to the vortex displacement the peak value of tw is more sensitive to the changing of RePL for shear-thickening rather than for shearthinning fluids. For RePLZ50 the maximum of tw for nZ1.5
678
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
(i) 1.0
(i) 1.0
Re=50
0.8
0.8
Re=100
n=3
Re=200
n=2
Re=300 0.6
n=1.5
0.6
Re=500
y
y
n=1 n=0.75
0.4
0.4
0.2
0.2
n=0.5
0.0
0.0 –0.4
–0.2
0.0
0.2
0.4
0.6
0.8
-0.4
1.0
-0.2
0.0
0.2
0.4
0.6
0.8
1.0
u
u (ii) 0.4
(ii) 0.4 Re=50
n=3
Re=100
n=2 n=1.5
Re=200
0.2
Re=300
n=1
0.2
n=0.75
Re=500
n=0.5
v
v
0.0 0.0
–0.2
-0.2 –0.4
–0.6
-0.4 0.0
0.2
0.4
0.6
0.8
1.0
x
0.0
0.2
0.4
0.6
0.8
1.0
x
Fig. 18. Velocity profiles for Power-Law fluid flow for nZ1.5 and different values of RePL number. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
Fig. 19. Velocity profiles for Power-Law fluid flow at RePLZ100 and different values of n. (i) u-velocity profile along the vertical centreline of the cavity; (ii) v-velocity profile along the horizontal centreline of the cavity.
is markedly greater than for nZ0.5 but it ends up smaller for RePLZ500. It can also be seen that for high values of RePL the higher values of the peak of tw occurs neither for high nor for low values of n but for values close to nZ1.
arrangement of variables. The convection terms are discretised using the QUICK scheme in order to avoid the artificial diffusion caused by interpolation schemes of low order of accuracy. Non-Newtonian flow computations were performed using four different non-Newtonian models namely the Power-Law, Bingham, Casson and Quemada models. The code is validated against both Newtonian and non-Newtonian lid-driven cavity flow cases. Results from the current study were in very good agreement with results by previous researchers for both Newtonian and Power-Law fluid flow cases. From the parametric study regarding the SCP of the Bingham, Casson and Quemada models
5. Summary A numerical code for the computation of viscous nonNewtonian flows has been presented. The code is based on the pressure-correction method and a collocated
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
679
Table 3 Maximum of the tw distribution at the bottom of the cavity for flow of Power-Law model with nZ1.5, 1 and 0.5 and different values of RePL
(i) 1.0
0.8
0.6
y
Re=50 Re=100 Re=200 Re=300 Re=500
0.4
0.2
n=0.5 n=1 n=1.5
0.0 0.0
0.2
0.4
0.6
0.8
1.0
RePL
nZ1.5
nZ1
nZ0.5
50 100 200 300 500
0.0201 0.0102 0.0053 0.0037 0.0025
0.0142 0.0076 0.0048 0.0042 0.0041
0.0090 0.0051 0.0039 0.0037 0.0035
fluids the increasing of Re number causes higher velocity gradients and subsequent moving of the centre of the main vortex towards the lid; (ii) the more shear-thinning the fluid is, the closer to the upper right edge the centre of the main vortex takes place; (iii) the peak value of the dimensionless wall shear stress ðtw Þ at the lower wall of the cavity is more sensitive to the changing of Re for shear-thickening fluids and can be higher or lower than the corresponding value for shear-thinning fluids depending on the value of Re.
x (ii) 1.0
References 0.8
0.6
y
n=3 n=2 n=1.5 n=1 n=0.75 n=0.5
0.4
0.2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
x Fig. 20. Position of centre of main vortex for flow of Power Law model: (i) with nZ1.5, 1 and 0.5 and different values of RePL; (ii) at RePLZ100 and different values of n.
the following conclusions are drawn: (i) the increasing of the SCP is followed by shortening of the distance between the centre of the main vortex and the lid whereas the profile for both the u- and v-velocities is flattened; (ii) the behaviour of the Bingham model is less sensitive to the changing of its SCP than that of the Casson and Quemada models; (iii) the Casson and Quemada models exhibit very similar trends as to the changing of their SCP. From the parametric study regarding the shear-thickening- and shearthinning effects on the flow the following conclusions are drawn: (i) for both shear-thinning and shear-thickening
[1] Anderson JD. Computational fluid dynamics. Singapore: McGrawHill; 1995. [2] Bell BC, Surana KS. p-Version least squares finite element formulation for two-dimensional, incompressible, non-Newtonian, isothermal and non-isothermal flow. Int J Num Meth Fluids 1994;18: 127–62. [3] Botella O, Peyret R. Benchmark spectral solutions on the lid-driven cavity flow. Comput Fluids 1998;27:421–33. [4] Brasseur E, Fyrillas MM, Georgiou GC, Crochet MJ. The timedependent extrudate-swell problem of an Olroyd-B fluid with slip along the wall. J Rheol 1998;42:549–66. [5] Dailey LD, Pletcher RH. Evaluation of multigrid acceleration for preconditioned time-accurate Navier–Stokes algorithms. Comput Fluids 1996;25:791–811. [6] Ferziger JH, Peric M. Computational methods for fluid dynamics. Berlin: Springer; 1999. [7] Ghia U, Ghia KN, Shin CT. High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method. J Comput Phys 1982;48:387–411. [8] Grillet AM, Yang B, Khomami B, Shaqfeh ESG. Modelling of viscoelastic lid driven cavity flow using finite element simulations. J Non-Newt Fluid Mech 1999;88:99–131. [9] Leonard BP. A stable and accurate convective modelling procedure based on quadradic upstream interpolation. Comput Meth Appl Mech Engng 1979;19:59–98. [10] Mitsoulis E. Three dimensional non-Newtonian computations of extrudate swell with the finite element method. Comput Meth Appl Mech Engng 1999;180:333–44. [11] Mompean G. On predicting abrupt contraction flows with differential and algebraic viscoelastic models. Comput Fluids 2002;31:935–56. [12] Neofytou P, Drikakis D. Effects of blood models on flows through a stenosis. Int J Num Meth Fluids 2003;43:597–635. [13] Neofytou P, Drikakis D. Non-Newtonian flow instability in a channel with a sudden expansion. J Non-Newt Fluid Mech 2003; 111:127–50. [14] Oliveira PJ, Pinho FT, Pinto GA. Numerical simulation of non-linear elastic flows with a general collocated finite volume method. J NonNewt Fluid Mech 1998;79:1–43.
680
P. Neofytou / Advances in Engineering Software 36 (2005) 664–680
[15] Papanastasiou TC. Flow of materials with yield. J Rheology 1987;31: 385–404. [16] Patankar SV, Spalding DB. A calculation procedure for heat, mass and momentum transfer in three dimensional parabolic flows. Int J Heat Mass Transfer 1972;15:1787–806. [17] Pham TV, Mitsoulis E. Entry and Exit Flows of Casson Fluids. Can J Chem Engng 1994;72:1080–4. [18] Quemada D. Rheology of Concentrated Disperse systems III. General features of the proposed non-newtonian model. Comparison with experimental data. Rheologica Acta 1978;17:643–53. [19] Rahman MM, Miettinen A, Siikonen T, Modified SIMPLE. formulation on a collocated grid with an assessment of the simplified QUICK scheme. Numer Heat Transf, Part B 1996;30: 291–314. [20] Rhie CM, Chow WL. A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation. AIAA J 1983;21: 1525–32.
[21] Saramito P, Roquet N. An adaptive finite element method for viscoplastic fluid flows in pipes. Comput Meth Appl Mech Engng 2001;190:5391–412. [22] Smyrnaios DN, Tsamopoulos JA. Squeeze flow of Bingham plastics. J Non-Newt Fluid Mech 2001;100:165–90. [23] Stone HL. Iterative solution of implicit approximations of multidimensional partial differential equations. SIAM J Num Analysis 1968;5:530–58. [24] Syrja¨la¨ S. Finite-element analysis of fully developed laminar flow of Power-Law non-Newtonian fluid in a rectangular duct. Int Comm Heat Mass Transf 1995;22:549–57. [25] Tannehill JC, Anderson DA, Pletcher RH. Computational fluid mechanics and heat transfer. Washington, DC: Taylor and Francis; 1997. [26] Xue SC, Phan-Tien N, Tanner RI. Upwind with deferred correction (UPDC): an effective implementation of higher order convection schemes for implicit finite volume methods. J Non-Newt Fluid Mech 2002;108:1–24.